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Abstract 

We  study  the  propagation  of  waves  over  large  regions  of  space  with  smooth,  but 
not  necessarily  constant,  material  characteristics,  separated  into  sub-domains  by 
interfaces  of  arbitrary  shape.  We  consider  a  divide  and  conquer  approach  based  on 
wave  splitting  into  incoming  and  outgoing  waves.  We  assemble  the  overall  solution 
from  the  set  of  individual  solutions  to  an  auxiliary  problem  (AP).  The  AP  is  defined 
independently  for  each  sub-domain.  The  choice  of  the  AP  is  relatively  flexible;  it 
can  be  formulated  to  enable  an  easy  and  economical  numerical  solution.  Our  new 
method  uses  only  simple  structured  grids,  e.g.,  Cartesian  or  polar,  regardless  of  the 
shape  of  the  boundaries  or  interfaces.  In  the  regions  of  smoothness,  it  employs  high 
order  accurate  finite  difference  schemes  on  compact  stencils.  They  do  not  require  any 
additional  boundary  conditions  besides  those  needed  for  the  underlying  differential 
equation  itself.  Interfaces  not  aligned  with  the  grid  handled  by  Calderon’s  operators 
and  a  method  of  difference  potentials  [  ]. 

The  operator  of  Calderon  and  the  method  of  difference  potentials  have  a  number 
of  important  advantages;  it  easily  handles  curvilinear  boundaries,  variable  coeffi¬ 
cients  and  general  boundary  conditions  while  the  complexity  is  that  of  a  finite- 
difference  scheme  on  a  regular  structured  grid.  A  main  advantage  is  that  this 
methodology  provides  high  order  accuracy  and  overcomes  the  difficulties  inherent 
in  more  traditional  approaches. 
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Chapter  1 


Introduction 


We  consider  problems  that  involve  the  propagation  of  acoustic  or  electromagnetic 
waves  over  large  regions  of  space  with  smooth,  but  not  necessarily  constant,  material 
characteristics,  separated  by  interfaces  of  arbitrary  shape.  The  external  boundaries 
can  also  be  arbitrarily  shaped.  These  problems  play  a  central  role  in  imaging  (med¬ 
ical  and  other  types),  nondestructive  evaluation,  land  mine  detection,  active  control 
of  sound,  etc.  We  describe  the  system  governed  by  the  Maxwell  equations  coupled 
with  appropriate  initial  and  boundary  conditions.  The  acoustic  case  results  in  a 
similar  situation.  The  Maxwell  equations  are  given  by: 


Vxl 


0  (Faraday  s  Law), 


df) 

dt 


Vxl 


-7 


(Ampere' s  Law), 
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CHAPTER  1.  INTRODUCTION 


coupled  with  Gauss’s  law  V  •  B  =  0,  V  •  if  =  p,  where  =  trlf  is  electric  current 
density,  a  is  electrical  conductivity,  and  p  is  electric  charge  density.  For  linear 
materials  we  relate  the  magnetic  flux  density  vector  if  to  the  magnetic  field  vector 
if  and  the  electric  flux  density  vector  if  to  the  electric  field  vector  if  using  if  = 
pt,  D  =  elf,  where  e  is  the  dielectric  permittivity  that  describes  the  particular 
media  and  p  is  the  magnetic  permeability.  We  consider  discontinuities  in  the  media, 
i.e.  in  e. 

In  two  dimensions  the  set  of  6  equations  decouples  into  two  independent  sets  of 
3  equations  denoted  as  TM,  and  TE^  (transverse  magnetic  and  electric  fields). 


Consider  the  two  dimensional  TM2  mode  Maxwell  equations  in  a  lossless  mate¬ 
rial,  i.e.  a  =  0: 


dHx 

1  dEz 

dt 

y  dy 

dHy 

_  1 

dEz 

dt 

A4 

dx  ’ 

dE, 

_ 1 

(dHy 

dt 

£ 

\  dx 

dHx 
dy 

Differentiating  the  first  equation  in  y,  the  second  in  x  and  the  last  in  t,  assuming 
£,  p  are  time  independent,  one  gets 


d2Hx  =  d  ( 1  dEz  \ 
dtdy  dy  yy  dy  J  ’ 

<  92 Hv  =  9_  (ldEz\ 
dtdx  dx  yy  dx  J  ’ 

d2Ez  _  1  ( d2Hy  d2Hx  \ 

dt 2  £  y  dxdt  dydt  J  ' 


Substituting  the  first  two  equations  into  the  third  one,  we  get  the  wave  equation 
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for  Ez: 


d2Ez  _  1  1  0EZ 


dt 2  e  \dx  \fx  dx 


Since  for  most  materials  y  is  constant  while  e  is  not,  we  rewrite  the  last  equation  as 


where  c2(x,y)  =  £,  ■  When  we  apply  the  Fourier  transform  in  time  the  wave 


equation  becomes  the  Helmholtz  equation  A u  +  k2u  =  0  where  k 2 


is  the 


wavenumber,  k  =  k(x,y). 

In  physical  applications,  the  quantity  c(x,  y)  may  be  discontinuous  and  therefore 
the  wavenumber  may  be  only  a  piecewise  continuous  function.  The  propagation  of 
waves  across  media  with  material  discontinuities  is  encountered  in  a  wide  variety  of 
settings.  For  example,  classical  problems  of  electromagnetic  scattering/transmission 
often  involve  sharp  variations  of  material  properties.  These  problems  appear  in 
applications  that  range  from  radar  imaging  to  telecommunication  devices. 

Our  discussion  here  concerns  the  numerical  solution  of  the  Helmholtz  equation 
for  domains  where  the  boundaries  and  interfaces  do  not  necessarily  conform  to 
the  discretization  mesh.  The  material  properties  are  assumed  smooth  between  the 
interfaces,  whereas  at  the  interfaces  they  may  undergo  jumps.  In  Section  1.1  we 
describe  the  general  problem  and  introduce  the  Helmholtz  equation  in  different 
coordinate  systems.  In  Section  1.2  we  review  existing  numerical  methods  while 
in  Chapter  2  we  explain  our  approach  in  solving  problems  with  discontinuities. 
In  Chapter  3  we  present  details  of  our  implementation  with  representative  results 
followed  by  our  conclusions  in  Chapter  4. 
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1.1  Helmholtz  equation 

The  Helmholtz  equation,  named  for  Hermann  von  Helmholtz,  is  the  elliptic  partial 
differential  equation 

A  u  +  k2u  =  F,  (1.1) 

where  A  =  V2  is  the  Laplacian,  u  =  u(x)  is  the  scalar  unknown  field,  e.g.,  acoustic 

pressure  or  linearly  polarized  electric  field  (a;  E  Mn),  and  F  =  F(x)  is  the  source 

term,  which,  if  present,  will  always  be  compactly  supported  in  Mn.  The  quantity 

2 

k  =  k(x)  in  (1.1)  is  the  wavenumber,  k 2  =  k^u2  and  k%  =  where  ojq  is  the  fixed 
carrier  frequency,  c  is  the  propagation  speed  in  the  unobstructed  medium  (also  fixed, 
as  the  speed  of  sound  in  ambient  fluid  at  a  constant  temperature  or  the  speed  of 
light  in  vacuum),  and  v  =  v(x)  is  the  refraction  index.  The  physical  interpretation 
of  v{x)  is  the  ratio  of  the  reference  speed  c  to  the  actual  propagation  speed  at  a 
given  x.  The  function  v{x),  and  hence  k(x),  can  be  only  piecewise  continuous. 
The  Helmholtz  equation  is  used  to  model  a  variety  of  important  physical  processes 
and  phenomena  in  acoustics  and  electromagnetism.  In  this  thesis  we  consider  two 
dimensional  problems  with  material  discontinuities  and  boundaries  that  are  not 
aligned  to  the  numerical  grid. 

Consider  an  incident  wave,  n^mc\  impinging  on  an  arbitrary  body  fii,  see  Figure 
1.1.  It  generates  a  transmitted  wave,  tt(trans);  and  partially  gets  reflected,  ■u(scat). 
It  could  be  an  ultrasound  wave  scanning  a  kidney,  an  embryo,  a  carrier  wave  of  a 
cellular  phone,  a  Wi-Fi  radio  signal  passing  through  a  wall  or  the  head  of  a  human. 
In  the  frequency  domain,  such  a  scenario  is  described  using  the  Helmholtz  operator 
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Lq  =  A  +  kq  in  the  domain  Llq  where  q  6  0, 1  as  follows 


(  LqU  = 

0 

x  e  q0, 

(1.2a) 

{  Liu  = 

F(x) 

x  6  fii, 

(1.2b) 

F(x)  is  the  source  term,  u(x)  =  n^mc^(a:) 

_l_  ^Iscat)^ j  for  each  x  e  Hq- 

u(inc)  = 

^—ik(x cos  d+y sm  o)  jg  a  gjven  piane  wave  impinging  from  an  incident  angle  9.  Thus, 
problem  (1.2)  is  driven  by  incident  wave  tt(inc)  and  the  source  term  F(x).  Across 
the  interface  T  between  the  two  subdomains,  Qi  and  Oo,  one  typically  requires  that 
the  function  and  it’s  first  normal  derivative  be  continuous.  Such  problems  may  have 
multiple  solutions.  To  ensure  uniqueness  one  additionally  requires  that  the  scattered 
field  ?Ascat)  has  no  incoming  components.  This  is  guaranteed  by  the  (n-dimentional) 
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Sommerfeld  radiation  condition: 


du(scat\x)  fscat).  .  /. 

- ;  +  ik0u^cat>(x)  =  o 

a \x\  \ 


l  —  n 

X\  2 


as  \x\ 


oo. 


(1.3) 


Problem  (1.2)  can  be  generalized  for  several  bodies  {0(?}^=1  C  Ho  which  are 
either  mutually  disjoint  or  share  an  interface.  The  problem  is  then  defined  as  fol¬ 
lowing: 

L0u  =  0  x  G  Ho, 

L\u  =  Fi(x)  x  G  Hi, 

(1-4) 


Tn'u  —  En (x)  x  G  Hn, 


driven  by  tt(inc)  and  a  set  Fj,  i  =  1,  2, . . . ,  n.  Across  every  interface  one  may  require 
that  the  function  and  it’s  first  normal  derivative  are  continuous.  One  can  also  define 
other  interface  conditions,  those  still  can  be  handled  by  the  method  described  below, 
see  equation  (2.18).  At  infinity  the  Sommerfeld  condition  is  required  for  uniqueness. 


On  the  other  hand,  when  one  solves  only  an  exterior  (Ho)  or  an  interior  (Hi) 
problem  then  only  one  boundary  condition  should  be  provided.  For  an  exterior 
problem,  in  addition  to  the  Sommerfeld  condition,  one  often  requires  that  u^scat^  = 
_M(mc)  or  riiscat)  =  — Un nc)  on  T  which  is  the  Dirichlet  u|r  =  0  and  Neumann 
Un\r  =  0  scattering  problem,  respectively.  un,  !iiscat*  and  are  the  first  normal 
derivatives  of  u ,  u^sca,t\  and  t/inc)  respectively. 


Mathematical  problems  are  often  described  better  in  some  specific  coordinate 
system  while  in  another  system  of  coordinates  the  problem  may  have  a  less  conve- 


32 


1.1.  HELMHOLTZ  EQUATION  29 

nient  representation.  In  boundary  value  problems  one  may  change  the  coordinate 
system  to  match  the  shape  of  the  boundary,  e.g.  use  polar  coordinates  for  circular 
domains  or  elliptical  coordinates  domains  shaped  as  ellipses.  However,  when  the 
shape  becomes  more  complicated  the  change  of  coordinates  may  become  a  major 
difficulty.  Furthermore,  choosing  a  complicated  coordinate  system  may  lead  to  an 
poor  quality  grid  or  singular  points,  thus  ruining  the  advantage  of  changing  coordi¬ 
nates.  Hence,  there  is  a  tradeoff  between  the  complexity  of  the  grid  and  complexity 
of  the  problem. 

For  a  generally  shaped  body,  constructing  a  body  fitted  coordinate  system  needs 
to  be  done  numerically.  Achieving  this  with  high  accuracy  may  be  difficult.  For  a 
complicated  body,  especially  in  3D,  it  may  be  impossible  to  construct  a  simple  body 
fitted  coordinate  system.  A  partial  remedy  may  be  to  use  the  multi-block  overlap¬ 
ping  grids,  also  known  as  chimera  grids,  see,  e.g.,  [  ].  Those  grids,  however,  do 

not  completely  remove  the  difficulty  associated  with  fitting  the  grid  to  a  curvilinear 
boundary.  They  rather  partially  alleviate  it,  because  the  grid  no  longer  has  to  be 
fitted  to  the  entire  boundary.  Instead,  different  blocks  serve  different  fragments  of 
the  boundary,  and  instead  of  point  matching,  the  grids  are  allowed  to  overlap  on 
some  common  regions,  which  simplifies  their  generation.  Yet  another  alternative  is 
to  use  the  finite  element  method  based  on  an  unstructured  grid,  see  Section  1.2. 

In  this  work  we  use  a  finite  difference  scheme  with  a  regular  grid  regardless  of 
the  shape  of  the  body,  and  still  obtain  high  order  accuracy  of  the  numerical  solution. 
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1.1.1  Helmholtz  equation  in  different  coordinate  systems 


We  next  recast  the  Helmholtz  equation  two  non-Cartesian  coordinate  systems.  The 
main  difficulty  is  the  Laplacian,  which  is  defined  as  divergence  of  the  vector  gradient, 
i.e.  Art  =  div(gradu)  ;  one  denotes  div  =  V-  and  grad  =  V  to  get  A  =  V  •  V  =  V2. 
In  the  Cartesian  coordinates  one  gets 


A  =  V- 


d2 


Things  are  more  complicated  in  polar  coordinates,  which  are  given  as 


x  =  r  cos  9, 

< 

y  =  r  sin#, 


where  r  >  0  is  a  radius  and  0  <  9  <  2ir  is  the  angular  coordinate.  Using  the  chain 
rule  one  gets 
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Such  a  complicated  derivation  is  usually  not  convenient.  Fortunately,  bet¬ 
ter  methods  exist.  Consider  an  orthogonal  curvilinear  coordinate  system  y  = 
y(x),x  G  Mn.  Define  scale  factors  (also  called  metrics  or  Lame  coefficients)  as 
h2j  =  £m= i  (f^)  •  Then  the  vector  gradient  becomes  V  =  (^g^,  ■  ■  • , 
and  the  Laplacian  becomes 

^  _  1  d  ( hp  d  \ 

n'/=i  hj  ^  dxm  V  hm  dxm)  ' 

Returning  to  polar  coordinates  one  gets 

hr 

hg 
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and  the  Laplacian  is: 


A  = 


1 

hr  hg 
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a 

8r 


h^d_ 

hr  dr 
d2 
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dd 


d  d2  1  d2 

dr  +  dr2  r  dd2 


hr  d 
hg  dd 
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1  d 


d_ 

dr 

1  d2 
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1  d_ 
r  dd 


^  r  dr  +  r2  dd2 


ip  =  7r/2 


Figure  1.2:  Elliptical  Coordinates 

Elliptical  coordinates  are  given  by  x  +  iy  =  d  cosh (77  +  icp)  or,  equivalently, 


x  =  d  cosh  r]  cos  <p 

< 

y  =  d  sinh  r/  sin  tp 


x  =  a  cos  tp 
y  =  b  sin  tp 


(1.5) 


where  d  is  semi-focal  distance,  rj  is  an  elliptical  radius,  so  that  each  value  of  rj  defines 
an  unique  ellipse  with  d  =  \/a2  —  b2,  a  =  d  cosh  77,  and  b  =  d sinh  r/  are  the  major 
and  minor  semi-axes,  respectively,  and  ^  =  cos2  p  +  sin2  tp  =  1. 
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In  this  case  we  have 


and  consequently, 


1.2  Numerical  approximation  of  differential  equations 

Finite  difference  (FD)  methods  were  historically  the  first  methodology  for  the  nu¬ 
merical  solution  of  differential  equations  [  >4,  4!  ].  They  still  remain  a  very  powerful 
tool,  and  for  smooth  solutions  on  regular  grids  lead  to  inexpensive  and  efficient  algo¬ 
rithms.  Their  primary  disadvantage  is  in  dealing  with  more  complicated  geometries 
and  solutions  with  low  regularity.  In  particular,  when  the  boundary  is  not  aligned 
with  the  grid  staircase  meshing  doesn’t  provide  the  required  accuracy.  For  example, 
consider  a  circular  body  Pi  in  Cartesian  coordinates,  see  Figure  1.3.  Denote  the 
circular  boundary  shape  T  =  90]  and  let  (x,  y )  G  T.  Let  vl,j  be  the  approximate 
value  of  the  solution  v  at  nodes  Xi  =  x  +  5x  and  yj  =  y  +  Sy,  i.e.  Vij  ~  v(xi,  yj). 
Let  b(x,y)  be  the  boundary  value  given  on  curve  T  and  assume  ( Xi,yj )  ^  T.  The 
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Figure  1.3:  Circular  body  in  Cartesian  coordinates,  staircase-mesh  7. 
staircase  approach  states  Vij  =  b(x,y).  However,  since 

v(xi,  yj)  =  b(x,  y )  +  Vv(x,  y )  •  (5x,  by)  +  0(bx2  +  by2), 


and  Vu  is  not  available,  the  boundary  data  is  approximated  by  a  first  order  method. 
For  additional  discussion  about  staircase  meshing  and  it’s  disadvantages  see  [  ]. 

The  immersed  boundary  method  (IBM)  of  Peskin  [  ]  requires  a  modification  of 

the  governing  equations  to  treat  the  geometric  irregularities,  and  a  smoothed  ap¬ 
proximation  of  the  5  function  to  treat  the  discontinuity  at  the  interface,  to  achieve 
first  order  accuracy.  It  has  been  improved  upon  with  the  immersed  interface  method 
(IIM)  of  LeVeque  and  Li  [  ]  (also  see  their  book  [38]).  They  enforce  the  jump  at  the 

interface  directly  in  the  finite  difference  scheme.  The  implementation  suffers  from 
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increased  complexity,  but  second  order  accuracy  can  be  achieved.  Later  Zhang  and 
LeVeque  [7  ]  used  fictitious  points  to  modify  the  discrete  linear  system  to  account 
for  the  correct  position  of  the  interface  and  the  proper  physical  interface  condi¬ 
tions.  Recently  Xu  [  ]  developed  efficient  and  stable  boundary  condition  capturing 

immersed  interface  method  for  simulating  a  flow  around  an  object.  The  method 
have  almost  a  second  order  accuracy  for  the  velocity  and  above  first  order  for  the 
pressure.  Johansen  and  Colella  [  ]  used  Embeded  boundary  to  solve  Poison  equa¬ 

tion  with  variable  coefficients  and  Dirichlet  boundary  condition  on  irregular  domain 
using  Cartesian  grid  with  second  order  accuracy.  They  extended  the  solution  into 
a  fictitious  domain  and  the  resultant  linear  system  is  non-symmetric.  However  it 
is  compatible  with  a  multigrid  and  adaptive  mesh  techniques  which  should  improve 
the  complexity.  These  techniques  though  aren’t  immediately  extendable  for  the 
Helmholtz  equation.  Recently,  Crockett,  Colella  and  Graves  [1  ]  used  this  method 
for  Poison  and  heat  equations  with  discontinuous  coefficients  and  reached  second 
order.  To  the  best  of  our  knowledge,  there  are  no  reported  uses  in  the  literature  of 
those  methods  for  anything  but  simple  Dirichlet,  Neumann,  or  interface  conditions 
(continuity  of  the  solution  and  its  normal  flux),  changing  the  boundary  condition 
requires  major  changes  to  the  algorithm  [  ],  and  extension  to  higher  than  second 

order  accuracy  is  not  straightforward.  The  ghost-cell  immersed  boundary  method 
(GCIBM)  [  ]  employes  extrapolation  to  impose  the  boundary  condition  implic¬ 
itly;  and  the  immersed  boundary  projection  method  (IBPM)  [  ]  treats  the  no-slip 

boundary  conditions  along  with  the  immersed  boundary  as  a  Lagrange  multiplier. 
All  these  methods  are  difficult  to  generalize  to  high  order  accuracy. 


A  high  order  method  (up  to  sixteenth-order  accuracy)  for  elliptic  equations,  for 
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a  matched  interface  and  boundary  method  (MIB),  was  introduced  by  Zhou  et  al. 

74,  7  ].  The  high-order  interface  conditions  are  implemented  by  repeatedly  match¬ 
ing  the  interface  conditions  across  the  given  interface  using  low-order  numerical  in¬ 
tegration  rules.  A  special  procedure  is  proposed  to  determine  the  accurate  fictitious 
values  required  for  the  high-order  scheme.  This  method  is  meant  to  treat  interface 
curves  that  are  not  aligned  to  the  grid,  jumps  in  coefficients,  and  singular  sources. 
However,  all  the  examples  presented  in  [  1,  7  ]  reduce  to  a  problem  with  only  singu¬ 
lar  sources.  Six  new  versions  of  IIM  of  fourth  order  accuracy  were  provided  by  Zhong 
].  He  used  two  different  polynomials  on  both  sides  of  the  interface  and  enforced 
that  the  two  polynomials  satisfy  two  interface  conditions.  Zhong  used  a  wide  stencil 
that  requires  additional  purely  numerical  boundary  conditions.  The  only  examples 
provided  were  the  Poisson  equation  with  a  singular  source  or  equivalent. 

An  interesting  solver  was  presented  by  Abarbanel  and  Ditkowski  [  ].  They  con¬ 
sidered  a  diffusion  equation  in  one  and  two  dimensions  on  a  domain  with  a  body 
whose  boundary  points  do  not  coincide  with  the  nodes  of  a  rectangular  mesh.  In 
order  to  reach  the  fourth  order  accuracy,  energy  methods  were  used  in  conjunction 
with  simultaneous  approximation  terms  (SAT).  However,  the  resulting  linear  system 
should  be  negative  definite,  which  is  not  the  case  for  the  Helmholtz  equation. 

The  finite  element  method  (FEM)  [  8,  13,  7]  and  its  extensions  (GFEM,  XFEM, 
discontinuous  enrichment  [  2,  23]),  as  well  as  discontinuous  Galerkin  methods  [  ], 

are  also  well  established  and  powerful.  Their  strength  is  in  dealing  with  complex 
geometries  and  low  regularity  of  the  solutions. 


In  practical  problems  of  wave  propagation  though,  especially  in  3D,  both  FD 
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and  FEM  have  serious  limitations  because  of  their  relatively  high  “points-per- 
wavelength”  requirement,  as  well  as  numerical  dispersion  and,  more  generally,  nu¬ 
merical  pollution  [  5,  Section  4.6.1],  [5,  ].  The  numerical  phase  velocity  of  the 

wave  in  these  methods  depends  on  the  wavenumber  k,  so  a  propagating  packet  of 
waves  with  different  frequencies  gets  distorted  in  the  simulation.  Furthermore,  the 
numerical  error  depends  strongly  on  the  wavenumber  k,  see  [  15],  and  this  kind  of 
error  is  inherent  in  FEM/FD.  The  error  behaves  like  hPkp+1  where  p  is  the  order  of 
accuracy  of  the  scheme.  So  the  number  of  points  per  wavelength  needed  for  a  given 
accuracy  grows  like  k1^.  Hence,  for  higher  order  accurate  schemes  the  pollution 
effect  is  reduced. 

A  high  order  approximation  can  be  built  for  arbitrary  boundaries  using  FEM, 
but  only  in  fairly  sophisticated  and  costlier  algorithms  with  isoparametric  elements 
].  In  discontinuous  enrichment  /  discontinuous  Galerkin  methods  and  GFEM, 
high  order  accuracy  also  requires  additional  degrees  of  freedom,  which  entails  an  ad¬ 
ditional  computational  cost.  These  additional  degrees  of  freedom  lead  to  expanded 
approximating  spaces  which  are  capable  of  approximating  a  very  broad  class  of 
functions  that  may,  in  principle,  have  irregularities  anywhere  in  the  computational 
domain.  This  is  a  significant  advantage  in  problems  of  great  geometrical  and  phys¬ 
ical  complexity.  However,  for  simpler  problems  with  smooth  solutions  much  more 
targeted  and  economical  approximations,  in  narrower  functional  spaces,  are  available 
using  FD  methods. 

A  FD  approach,  on  the  other  hand,  does  not  introduce  additional  unknowns  per 
grid  node  and  thus  remains  inexpensive.  It,  however,  requires  a  higher  regularity 
of  the  solution  to  guarantee  consistency,  and  may  also  need  a  wider  stencil,  which 
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complicates  the  boundary  conditions  (more  precisely,  requires  additional  “purely 
numerical”  boundary  conditions). 

Schemes  known  as  Collatz  “Mehrstellen”  [  ],  equation-based  and  related  com¬ 

pact  schemes  [55,  66,  69,  33,  3,  ,  S,  ],  as  well  as  Trefftz-FLAME  [  ],  don’t  expand 
the  stencil  and  therefore  don’t  require  additional  nonphysical  boundary  conditions. 
Furthermore,  as  opposed  to  classical  FEM  which  expands  the  approximating  space, 
these  methods  reduce  it  to  the  class  of  solutions  rather  than  the  much  broader  class 
of  generic  sufficiently  smooth  functions.  This  does  not  imply  any  loss  of  generality, 
because  according  to  the  Lax  theorem,  for  convergence  one  does  not  need  to  have 
consistency  for  any  functions  except  the  solutions. 

1.2.1  Compact  Equation— Based  Schemes 

We  next  present  the  compact  scheme  used  hereafter  in  our  simulations.  Consider 
Uij  =  u(xi,  r/j),  then  the  second  order  accurate  approximation  to  the  second  deriva¬ 
tive  in  x  yields 

n  _  ui+i,j  ~  2Uij  +  Ui- ij 

-L'xxU/iJ  —  ^  2 

Using  Taylor  series  one  can  show  that 

h2 

D xx ^ i ■']  —  'U’xx  T  TT^im  T  O  (j~l  )  . 


(1.6) 
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To  obtain  a  4th  order  approximation  one  eliminates  uxxxx  using  the  two  dimensional 
Helmholtz  equation  —uxx  =  uyy  +  k2u  —  F.  Differentiating  this  equation  we  obtain 

uXXXx  =  uyyxx  T  (k  uj xx  Fxx  =  DyyDxxu  T  Dxx  (^k  u  F'j  -T  0(/i  ) 

where  Dyyuid  =  a‘^+1~2u^j+Ui^-1 .  Thus,  the  fourth  order  compact  scheme  is  given 
by 

h2 

DxxUij  —  DxxUij  —  ( Oyy Dxxv  T  Dxx  ( k  u  T) ) 

=  uxx  +  O  (/i4)  .  (1-7) 

Using  a  similar  approach  for  the  derivative  in  y,  one  gets  the  scheme  of  Singer  and 
Turkel  [55] 

^ DXX  +  Dyy  +  k~  +  —  (/l^  +  h y)  DXXDyy^j  U  +  ~  (jl^D  XX  +  flyDyy )  (A)2u) 

=  +  —  ( hxDxx  +  hyDyy^j  F.  (1-8) 

Compact  fourth  order  accurate  schemes  for  the  more  general  Helmholtz  -type 
equations  with  variable  coefficients  have  been  obtained  in  [8,  9].  A  sixth  order 
scheme  for  constant  coefficients  is  constructed  in  [  ].  A  similar  compact  sixth 

order  scheme  for  the  case  of  a  variable  wavenumber  k  =  k(x )  is  presented  in  [  ]. 

A  distinctive  feature  of  the  compact  equation-based  schemes  is  that  they  exploit 
two  stencils.  The  first  one  applies  to  the  left-hand  side  of  the  equation,  i.e.,  it 
operates  on  the  unknown  solution.  The  second  one  applies  to  the  right-hand  side  of 
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the  equation,  i.e.  ,  it  operates  on  the  given  data  (source  terms).  The  equation-based 
and  similar  high  order  schemes  reduce  pollution  while  keeping  the  treatment  of  the 
boundary  conditions  simple.  Since  the  order  of  the  resulting  difference  equation  is 
equal  to  the  order  of  the  differential  equation,  no  nonphysical  boundary  conditions 
are  required. 

1.2.2  Reduction  to  integral  equations 

In  traditional  boundary  element  methods  (BEM),  linear  boundary  value  problems 
are  transformed  into  integral  equations  with  respect  to  equivalent  boundary  sources, 
which  are  subsequently  discretized.  Practical  applications  of  such  methods  date 
back  to  the  1960s.  They  impose  no  limitations  on  the  shape  of  the  boundary  and 
automatically  account  for  the  correct  far  field  behavior.  There  are,  however,  several 
major  disadvantages: 

•  Full  matrices  —  in  contrast  with  the  sparse  FD  and  FEM  matrices.  (Cases  of 
quasi-sparse  integral  equations,  due  to  the  rapid  decay  of  Green’s  functions  in 
space,  are  exceptional,  [46]). 

•  A  relatively  narrow  treatment  of  the  boundary  conditions.  Care  must  be  exer¬ 
cised,  on  a  case-by-case  basis,  in  the  choice  of  the  equivalent  boundary  sources, 
so  that  the  resulting  Fredholm  equation  is  of  the  second  kind  (well-posed) 
rather  than  first.  Moreover,  mixed  (Dirichlet/Neumann,  etc.)  or  less  standard 
(Robin,  etc.)  conditions  require  a  special  development. 

•  Singular  integral  kernels  ( the  most  serious  disadvantage  in  practice).  Imme- 
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diately  at  the  boundary  points,  the  kernel  singularity  can  usually  be  handled 
analytically,  and  the  fields  remain  bounded  as  long  as  the  surfaces  are  smooth. 
However,  for  points  in  the  vicinity  of  a  surface,  the  evaluation  of  the  integral 
is  problematic,  as  analytical  expressions  are  usually  unavailable  and  numerical 
quadratures  require  extreme  care. 

•  Limitation  to  constant  coefficients,  i.e. ,  to  homogeneous  media,  —  these  meth¬ 
ods  require  explicit  knowledge  of  the  fundamental  solution  of  the  correspond¬ 
ing  differential  operator. 

Significant  progress  in  Fast  Multipole  Methods  (FMM)  [  5,  12,  ,  i,  ]  has 

helped  alleviate  the  first  disadvantage  of  boundary  methods.  FMM  accelerates  the 
computation  of  fields  due  to  distributed  sources  —  or  equivalently,  matrix-vector 
multiplications  for  the  dense  system  matrices.  But  the  second  and  especially  the 
third  disadvantages  are  more  difficult  to  overcome,  whereas  the  fourth  one  can  only 
be  addressed  by  switching  to  much  less  efficient  volume  integral  methods  or  finding 
the  fundamental  solution  by  numerical  means,  which  is  both  expensive  and  leads  to 


additional  errors. 
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Chapter  2 


Problems  with  non-aligned 
interfaces 


In  this  chapter  we  describe  our  main  methodological  approach  for  solving  problem 
(1.2).  We  consider  a  divide  and  conquer  approach  based  on  wave  splitting  into 
incoming  and  outgoing  waves.  We  assemble  the  overall  solution  from  the  set  of 
individual  solutions  to  the  auxiliary  problem  (AP).  The  AP  is  defined  independently 
for  each  sub-domain  Qq,  see  (1.2).  The  choice  of  the  AP  is  relatively  flexible;  it  can 
be  formulated  so  as  to  enable  an  easy  and  economical  numerical  solution.  Our  new 
method  uses  only  simple  structured  grids,  e.g.,  Cartesian  or  polar,  regardless  of 
the  shape  of  the  boundaries  or  interfaces.  In  the  regions  of  smoothness,  it  employs 
high  order  accurate  finite  difference  schemes  on  compact  stencils,  see  Section  1.2.1. 
They  do  not  require  any  additional  boundary  conditions  besides  those  needed  for 
the  underlying  differential  equation  itself. 
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We  first  illustrate  the  key  components  of  the  numerical  methodology  that  we 
propose  with  a  one  dimensional  example  and  then  generalize  this  simple  one  di¬ 
mensional  model  of  wave  splitting  into  a  multidimensional  one  using  Calderon’s 
projections,  see  [  ]. 

Consider  the  following  one-dimensional  problem 

uXx  +  klu  =  0  x  <  0, 

(2.1) 

uxx  +  kfu  =  0  x  >  0, 

driven  by  an  incoming  wave  u^llic^  =  Aelk°x,  which  propagates  through  the  domain. 
It  generates  the  transmitted  wave  u\  =  TelklX  for  x  >  0  and  partially  gets  reflected 
and  produces  uq  =  Re~lk°x  for  x  <  0,  see  Figure  2.1. 


Figure  2.1:  ID  transmission  problem. 
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To  find  the  amplitudes  T  and  R  one  requires  that  u(x)  and  u'(x)  be  continuous  at 
x  =  0.  This  yields  two  linear  algebraic  equations:  A  +  R  =  T  and  kpA  —  kpR  =  k\T , 
which  provide  the  Fresnel  reflection  and  transmission  coefficients: 


T 


=  A 


2fc  o 

k0  +  h ' 


R  =  A 


kp  -  h 
kp  +  ki  ’ 


(2.2) 


However,  one  can  solve  this  problem  without  introducing  unknown  amplitudes 
and  instead  employ  the  relation  that  describes  the  entire  family  of  waves:  right 
traveling  v!  —  ik\u  =  0  and  the  left  traveling  v!  +  ikpu  =  0.  In  particular,  the  trans¬ 
mitted  wave  u\  satisfies  the  right  traveling  equation.  Likewise  up  +  tthnc)  satisfies 
the  inhomogeneous  left  traveling  equation  on  x  <  0.  So  we  can  rewrite  (2.1)  as  the 
following  system  of  equations 


v!  —  ik\u  =  0 

v!  +  ikpu  =  2ikpAelk°x 


x  >  0, 
x  <  0, 


(2.3) 


together  with  the  requirement  that  u(x),u\x )  be  continuous,  which  yields: 


m(0)  =  A 


2kp 

kp  +  k\  ’ 


u'(0)  =  A 


2ik\kp 
kp  +  k\  ‘ 


(2.4) 


Since  T  =  A  +  R  =  u{ 0),  we  arrive  at  the  same  answer  (2.2). 

The  system  of  equations  (2.3)  shows  that  if  the  impinging  wave  changes  (new 
amplitude  A),  then  the  problem  does  not  have  to  be  solved  all  over  again,  because 
only  the  right-hand  side  of  inhomogeneous  part  changes.  In  one  dimension  the 
difference  is  negligible.  However,  in  multiple  space  dimensions  a  methodology  built 
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around  the  same  idea  proves  very  useful.  For  example,  in  inverse  problems  one  may 
cheaply  use  multiple  impinging  waves  or  different  incident  angles  to  get  additional 
information  to  improve  the  accuracy. 

In  Section  2.1  we  present  the  theoretical  part  for  Calderon’s  boundary  equations 
with  projections  which  replaces  the  system  of  equations  (2.3).  We  conclude  chap¬ 
ter  Chapter  2  describing  the  discrete  counterparts  of  the  boundary  equations  with 
projections  developed  by  V.S.  Ryaben’kii,  see  [50],  in  Section  2.3.  The  pseudocode 
of  the  main  algorithm  for  solving  (1.2)  is  described  in  Section  2.3.1,  followed  by 
examples  and  additional  explanation  in  Sections  2.3.2  -  2.3.7. 


2.1  Calderon’s  Potentials 

In  this  section  we  present  the  Calderon  potentials,  which  yield  the  wave  splitting 
described  above,  see  [  ].  All  the  incoming  and  outgoing  waves  for  a  given  Qq 

(analogues  of  the  right  and  left  traveling  waves  in  ID  example,  see  Figure  2.1)  belong 
to  the  image  (range)  and  kernel  (null  space),  respectively,  of  the  corresponding 
projection  operator.  We  consider  time-harmonic  wave  propagation  with  a  piecewise 
continuous  index  of  refraction. 

Let  Lq  denote  the  Helmholtz  operator  of  equation  (1.1)  with  k  =  kq,  where 
q  6  {0, 1},  and  the  geometry  as  depicted  in  Figure  1.1,  i.e.  we  consider  problem  (1.2). 
Let  Gq(x)  be  the  fundamental  solution  of  Lq.  Let  the  functions  £^( x )  and  £,n(x) 
belong  to  C2  ,  C1  respectively  for  *  e  T.  We  introduce  the  vector  =  (£7,£7/) 
which  we  call  the  vector  density.  A  generalized  potential  of  Calderon  type  with  the 
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vector  density  £r  is  defined  by 

Pnq£r{x)  =  crq j^{iI{y)d^-{x-y)-iII{y)Gq{x-y]sjdsy,  x  e  Dq,  (2.5) 

where  ^  denotes  the  normal  derivative  and  crq  =  (— l)9-1-1  changes  the  sign  so  that 
we  always  consider  a  normal  pointing  in  the  same  direction  regardless  of  the  domain 
Qq.  Given  the  solution  u(x)  to  the  problem  (1.2)  one  defines  a  vector  density 
up  =  (u,  f^)|r,  and  rewrites  Green’s  solution  to  the  problem  (1.2)  as: 


u(x) 


Pshur  +  Jn  Gi(x  -  y)F(y)dy  x  E  fii, 

< 

Pn0ur  x  e  n0. 


(2.6) 


Let  £r  =  (£7,£/7)  belong  to  the  space  S.  Inspired  by  the  one  dimensional 
example  we  wish  to  split  S  into  incoming  S+  and  outgoing  Si  waves  for  a  given  Qq 
such  that  their  direct  sum  is  S ,  i.e. 

S  =  S«©  Si.  (2.7) 

So  for  each  £  S  we  looking  for  E  S+  and  E  Si  such  that  +  £r  • 

Furthermore,  we  require  that  such  a  representation  be  unique.  One  implements  it 
by  a  projection  operator  whose  image  is  S+  and  it’s  kernel  is  Si.  We  stress  that 
while  the  space  S  doesn’t  change,  the  representation  (2.7)  changes  for  the  exterior  or 
interior  problem,  e.g.  S\_  ^  Si  in  general.  In  addition,  the  choice  of  the  projection 
operator  affects  the  representation  of  (2.7),  which  can  be  considered  as  changing 
the  projection  angle  onto  the  same  subspace  [  ]. 


51 


48  CHAPTER  2.  PROBLEMS  WITH  NON-ALIGNED  INTERFACES 

One  defines  such  a  projection  operator  using  Calderon’s  potential  (2.5).  This 
projection  is  known  as  the  Calderon’s  boundary  projection.  We  will  show  that  this 
projection  implies  a  wave  split,  see  discussion  in  Section  2.1.1.  We  first  define  the 
vector  trace  operator  Tr  ”=  (*>•  fair-  Thus,  Calderon’s  boundary  projection  is 
defined  as 

P&r=TrPnq£r.  (2-8) 

Note,  that  any  £p  satisfies  LqPn  =  0,  x  G  Llq,  and  therefore  Green’s  formula 
provides  Pn  £ r  =  Pfiq  Tr  Pnq£r  which  immediately  implies  that  Pp  is  a  projection 
since  (Pp)2  =  Pp. 

The  operator  Pp  has  an  important  property.  If  £p  defines  u{x)  =  Pnq€r  for 
which  Tr  u(x)  =  £p  ,  then  the  following  equation  holds 


P&r=£r.  (2-9) 

If  equation  (2.9)  holds,  then  u(x)  =  Pnq€r  satisfies  Lqu  =  0,  and  Tr  u  =  ^p. 
Finally,  we  have  proved  that  £p  satisfies  the  BEP  if  and  only  if  £p  =  Tru  for  which 
Lqu  =  0  [51,  53,  ].  We  call  equation  (2.9)  the  boundary  equation  with  projection 

(BEP). 

2.1.1  Wave  Split 

The  solutions  to  the  homogeneous  equation  Lqu  =  0  can  be  interpreted  as  incoming 
waves  for  the  domain  Llq,  because  they  have  no  (radiating)  sources  on  Clq.  In  order 
to  describe  how  Calderon’s  projection  splits  the  waves  [39],  we  first  recast  equation 
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(2.5)  as: 


PqA  r(*)  =  ‘ 


o, 


0 


fr  (£J(y)i£r(x  -y)~  £n(y)Gq(x  -  y)J  dsy ,  x  g  fiq, 


X  eRn\  fiq. 

(2.10) 


Next  we  define  Green’s  operator  on  Lq  as 


GqF{x)  =  < 


fnq  Gg(x  -  y)E(y)dy,  x  e  Qq, 

0  x  €  Mn  \  fig. 


(2.11) 


Thus,  one  expresses  the  exterior  solution  as  u°(x)  =  Pn0Up(x)  where 


uvr=  Tru°  =  u°, 


du°\ 
dn  ) 


and  the  interior  solution  as  the  ul{x)  +  G\  F(x)  =  Pq,  u^(x)  +  G\  F(x)  where 


uf  =  Tr  u 1  =  (  u1 , 


1  dul 


dn 


Hence,  the  solution  to  the  problem  (1.2)  is  given  by 


u(x)  =  u°(x)  +  v}(x)  +  G\  F(x),  x  G 


(2.12) 


We  next  denote  ur  =  uf  +  Up.  Note,  that  in  general  ur  does  not  satisfy  the 
BEP  unless  ur  =  Up.  This  means  that  the  operator  Pp£p  eliminates  that  part  of 
which  is  not  the  trace  of  the  solution  of  Lqu  =  0.  In  other  words  incoming  waves  to 
fig,  denoted  as  ,  belong  to  the  image  of  the  projection  Pp,  i.e.  £p  G  ImPp  while 
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outgoing  waves  to  Plq  satisfy  E  KerPp.  Therefore,  the  aforementioned  space  S 
satisfies  S  =  ImPp  ©  KerPp,  since  SJ+  =  ImPp  and  Si_  =  KerPp. 

The  operators  Pp  are  defined  in  such  a  way  so  that  if  there  is  no  discontinuity, 
i.e.  ko  =  k±  then  ImPp  =  KerPp  and  KerPp  =  ImPp,  and  so  there  are  no  reflections 
at  r.  However,  in  general  ImPp  ^  KerPp  and  KerPp  ^  ImPp,  which  means  that 
there  is  both  propagation  through  and  reflection  from  the  interface  T,  similar  to  the 
ID  case,  see  Figure  2.1. 

2.1.2  Reduction  of  the  problem  to  the  boundary 

For  an  exterior  domain  Do  the  solution  is  u  =  u+vync>  where  u  denotes  the  scattered 
field  satisfying  the  Sommerfeld  radiation  condition  (1.3).  Therefore,  due  to  the 
linearity  of  the  trace  operator,  the  density  for  the  exterior  domain  is  given  by 

£r  =  Tru  =  Tru  +  Tru ^  =kr+  £?nc)- 

To  satisfy  (1.3)  one  formulates  the  BEP  for  the  exterior  problem  for  the  scattered 
field: 

-fr£r  =  £r- 

However,  in  order  to  match  it  to  the  interior  part  one  complements  it,  for  the  total 
field,  by  adding  Pp  £p,1C')  +  £pmc^  to  both  sides  of  the  equation.  Hence,  assuming 
that  the  solution  and  its  first  normal  derivative  are  continuous  across  F,  we  can 
equivalently  reformulate  problem  (1.2)  as  the  following  system  of  BEP  [  ]  defined 

on  £r: 
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Pp£r  +  TrG\F  =£r, 

< 

P°tr  +  (I-P°)^nc)  =tr- 

System  (2.13)  can  be  thought  of  as  a  multi-dimensional  analogue 
Once  (2.13)  has  been  solved  for  £r,  the  solution  u{x)  is  given  by 


(2.13) 

of  system  (2.3). 


u(x) 


PnAr(x)  +  G\F(x)  x  E  fii, 

Pn „[fr  -^rmc)](®)  +  it{inc)(*)  x  e  ^o- 


(2.14) 


Note  that  the  application  of  the  trace  operator  Tr  reduces  system  (2.14)  back  to 
(2.13). 


2.1.3  Divide  and  Conquer 

One  of  the  important  generalizations  of  the  proposed  method  is  that  one  can  redefine 
(2.13)  in  such  a  way  so  that  it  is  possible  to  solve  it  as  two  separate  auxiliary  prob¬ 
lems  (AP)  one  for  the  interior  and  one  for  the  exterior  domain.  Such  an  approach 
has  several  advantages.  For  example,  one  can  solve  each  of  the  auxiliary  problems 
on  a  different  grid.  Another  example  is  that  one  can  change  the  Sommerfeld  approx¬ 
imation  without  recalculating  the  auxiliary  problem  for  the  interior  area.  However, 
the  main  advantage  is  the  relatively  simple  treatment  of  curvilinear  interfaces  T  as 
will  be  explained  in  the  next  section. 

Consider  an  arbitrary  function  w(x),  x  E  Rn,  that  satisfies  the  Sommerfeld 
radiation  condition  (1.3)  and  such  that  Trw  =  £r  =  (£7 ,£n)-  Generally  Lqw  /  0, 
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x  £  fig,  and  therefore  Green’s  formula  gives: 

w(x)  =  Pnq  Wr  +  Gg  Lq  W,  X  £  fig. 

Using  the  fact  that  Pfiq  wr  =  Pnq  £r(*),  we  get 

pnq€r(x)  =  w(x)  -  GgLgW,  xeflq.  (2-15) 

We  next  generalize  (2.15)  in  the  following  way:  let  Lq  be  defined  for  x  £  W 1 
(instead  of  a;  £  flq).  We  reformulate  Pnq  so  that  the  equation  (2.5)  becomes 

PnCI^r{x)  =  w(x)  -  Gq^{Lqw)\n^  ,  x  £  flq.  (2.16) 

The  potential  (2.16)  is  obtained  on  Mn  and  then  restricted  to  flq.  Furthermore,  it 
is  insensitive  to  the  choice  of  w  as  long  as  the  conditions  Trw  =  £r  and  (1.3)  hold. 
The  projection  Py  defined  through  (2.16)  is  the  same  as  one  defined  through  (2.5). 
The  importance  of  this  new  definition  is  that  it  does  not  contain  surface  integrals. 

This  further  generalization  is  the  key  feature  in  the  treatment  of  generally  shaped 
bodies  that  are  not  aligned  to  the  numerical  grid,  which  is  discussed  in  Section  2.3. 
Let  flq  be  a  regularly  shaped  expanded  domain  such  that  flq  C  flq.  Assume  that  Lq 
is  also  defined  on  flq,  and  that  Gq  is  the  corresponding  inverse  so  that  the  solution  u 
to  the  equation  Lqu  =  F  on  flq  is  given  by  u  =  GqF.  Hereafter,  we  assume  that  this 
solution  exists  for  any  given  F  defined  on  flq  and  is  unique  as  along  as  u  is  required 
to  satisfy  some  specially  chosen  boundary  conditions  at  dflq.  The  combination  of 
the  differential  equation  Lqu  =  F  and  these  boundary  conditions  will  be  referred  to 
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as  the  auxiliary  problem  (AP).  The  AP  can  be  constructed  so  that  it  is  easy  to  solve. 
In  particular,  the  evaluation  of  GqF  does  not  need  to  involve  any  singular  integrals. 
Instead,  the  AP  is  discretized,  and  its  solution  computed  using  finite  differences  on 
a  regular  grid  (see  Section  2.3).  Finally,  problem  (1.2)  is  solved  as  two  different 
BEPs  defined  on  £p: 

+  TrG\F  =#, 

(2.17) 

P^o+(/-pO)^mc)  =  £°, 

subject  to  the  condition 

A£  f  +  P£p  =  0,  (2.18) 

which  is  an  additional  generalization  to  (2.13)  where  A  =  B  =  1.  This  last  gener¬ 
alization  is  the  pathway  for  different  interface  conditions. 

We  denote  £p  =  (C;  •  £,!/)  where  q  E  {0, 1}  refers  to  the  interior  or  exterior  traces 
in  (2.17).  t^q1  are  the  solution  and  it’s  normal  derivative  on  the  interface  F, 
and  are  therefore  unknown.  In  order  to  solve  system  (2.17)  one  expands  both  £p, 
£p  in  some  basis,  e.g.  Fourier,  so  that  the  coefficients  of  the  expansion  become  the 
unknowns  of  the  problem.  Once  the  coefficients  are  known  and  £p,  £p  assembled, 
the  solution  is  computed  by  (2.14).  We  describe  this  approach  in  more  detail  and 
provide  examples  in  Section  2.3. 
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2.1.4  Boundary  conditions 

When  one  solves  only  the  exterior  or  interior  part  of  (2.17)  the  interface  condition 
on  T  becomes  a  boundary  condition  l-pu  =  4>.  One  recasts  it  as 


*r(^n1^f+GiF)  =  0  (2.19) 

for  an  interior  problem.  The  system  of  equations  (2.17),  (2.19)  is  then  solved  with 
respect  to  £p. 

For  an  exterior  problem  it  becomes 

k(Pn0tr  +  u(inc))  =  0  (2.20) 


and  the  system  of  equations  (2.17),  (2.20)  is  to  be  solved  with  respect  to  £r. 

The  operator  Zp  that  defines  the  boundary  condition  in  the  problem  can  be 
arbitrary,  ranging  from  very  simple  (e.g.,  Dirichlet  or  Neumann)  to  very  general 
(e.g.,  different  type  on  different  parts  of  T,  nonlocal,  etc.).  Systems  (2.17),  (2.19)  or 
(2.17),  (2.20)  are  still  equivalent  to  the  relevant  interior  or  exterior  problem. 

In  practice,  Calderon’s  potentials  and  projections  are  approximated  by  differ¬ 
ence  potentials  and  projections,  respectively,  see  [  ].  In  doing  so,  the  discrete 

counterparts  of  formula  (2.16)  are  developed  and  used,  so  that  one  never  needs 
to  approximate  singular  surface  integrals.  Instead,  one  needs  to  solve  a  discrete 
auxiliary  problem  that  can  be  chosen  convenient  for  a  numerical  solution. 


58 


2.2.  WELL-POSEDNESS  55 

2.2  Well-posedness 

We  assume  that  the  original  problem  (1.2)  is  well-posed,  i.e. ,  that  its  solution  exists, 
is  unique,  and  depends  continuously  on  the  data  F(x),  4>,  in  the  sense  of  appropri¬ 
ately  chosen  norms.  Then,  the  equivalent  problem  on  the  interface  (2.17),  (2.18) 
or  on  the  boundary  (2.17),  (2.19)  or  (2.17),  (2.20)  is  also  well-posed.  This  means 
that  if  the  BEP  is  perturbed,  then  the  solution  of  the  boundary  system  will  also  get 
perturbed,  and  the  perturbation  of  the  solution  will  be  bounded  in  the  appropriate 
norm  by  the  perturbation  introduced  into  the  BEP. 

For  instance,  consider  the  interior  homogeneous  case: 

Lu  =  0  on  11  and  l^u  =  (j)  on  T, 

for  which  the  equivalent  boundary  formulation  is 

-Pr£r-£r  =  0  and  lr{Pn€r)  =  4>- 
If  the  original  problem  is  well-posed,  then  ||u||  ^  c\\(/>\\,  and  consequently, 

IKrKdH# 

This  result  has  nothing  to  do  with  Calderon’s  operators  per  se,  it  holds  simply 
because  £r  =  Tru.  Let  ipr  be  a  perturbation,  so  that  instead  of  the  true  unperturbed 
boundary  problem  we  are  solving 


-Pr£r  —  £r  =  V’r  and  lr(Pn€r)  =  4>- 
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Then,  we  have 

where  the  constant  C  depends  on  |Pq||  and  ||Pr||,  but  does  not  depend  on  either 
4>  or  xpr-  The  proof  can  be  found  in  [  2,  Part  II,  Chapter  1],  It  is  based  on  splitting 
the  entire  space  of  traces  on  P  into  the  direct  sum:  ImPp  ©  KerPp  as  discussed 
above. 

2.3  Discrete  Calderon’s  Potentials 

We  now  develop  the  discrete  counterpart  of  the  theory  introduced  in  the  previous 
section.  To  make  the  presentation  easier  we  first  describe  it  in  relatively  simple 
algorithmic  terms  in  Section  2.3.1  followed  by  several  representative  problems  and 
their  solution  supported  by  BEP  theory  in  mathematical  and  algorithmic  terms  in 
following  sections.  We  stress,  despite  the  fact  that  we  only  provide  examples  with 
simple  regular  shapes,  the  theory  remains  unchanged  for  a  generally  shaped  body. 

2.3.1  Algorithm 

In  this  section  we  present  a  sequence  of  simple  algorithmic  steps  that  implements 
the  solution  of  problem  (1.2). 

Thus,  we  provide  pseudocode  of  the  main  procedures  required  to  implement  our 
method  and  give  several  algorithms  for  solving  the  problem  in  a  domain  not  aligned 
to  the  grid  including  the  transmission-reflection  algorithm  in  heterogeneous  media 
with  a  jump  at  the  interface  and  an  efficient  algorithm  for  multiple  impinging  waves. 
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The  main  advantage  of  our  method  is  that  for  each  sub-domain  Qq,  see  (1.2),  we 
compute  a  relatively  small  set  of  individual  solutions  to  a  simple  auxiliary  problem 
on  a  regular  domain  formulated  to  be  numerically  effective.  Individual  solutions  are 
used  to  assemble  the  solution  of  the  original  problem. 

More  precisely,  we  solve  the  system  of  BEPs  (2.17),  (2.18)  as  follows.  For  each 
BEP,  exterior  (q  =  0)  or  interior  ( q  =  1)  we  consider  an  expansion  of  =  (Cq^q1) 
in  some  basis  on  T,  e.g.  Fourier: 

Cqis)  =^2cnqbn{s),  (2.21a) 

and 

£,c/(s)  =  cn/,9Ms),  (2.21b) 

so  that  the  problem  is  reduced  to  finding  the  coefficients  chq ,  cj'q  that  satisfy  the 
condition 

A  (2  Cn°bn(s),  Cn’°bn(s))  +  B  (j]  Cn’1&n(s)>  ^Ms))  =  0  (2.22) 

obtained  by  substituting  expressions  (2.21)  into  formula  (2.18).  Each  individual 
solution  corresponds  to  a  basis  function  chosen  on  P  and  extended  to  it’s  grid  rep¬ 
resentation  7  by  means  of  the  operator  T  to  be  defined  in  Section  2. 3. 1.2. 

In  the  following  sections  we  define  7,  then  describe  the  extension  of  basis  func¬ 
tion,  and  then  describe  and  provide  pseudocode  of  the  discrete  counterpart  of 
Calderon’s  potentials  theory  followed  by  examples  of  solutions  to  several  problems. 
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2. 3. 1.1  The  grid  representation  of  the  boundary  shape 


Figure  2.2:  9-point  stencil  centered  at  node  Pij 


We  define  now  the  grid  representation  of  the  curve  T,  which  we  denote  7,  in  a 
following  way.  Let  12  be  a  smooth  body  in  the  domain  represented  by  regular  grid 
N,  e.g.  Cartesian.  Since  Llq  is  not  necessary  aligned  to  the  grid,  the  intersections 
between  N  and  T  may  be  an  empty  set,  therefore  we  define  7  to  be  a  set  of  grid 
points  surrounding  T.  Formally,  let  be  the  stencil  centered  at  node  pij,  e.g. 
the  9-point  stencil  on  Figure  2.2  and  let  N+,N“  be  non-empty  sets  defined  as 


n+=  U  NP«.  N"=  U 

where  q  E  {0, 1}.  Then  the  grid  representation  of  the  curve  I  ,  7,  is  given  as 


7  =  N+nr. 


(2.23) 


An  example  of  7  for  general  body  in  Cartesian  coordinates  is  shown  in  Figure  2.3. 
In  Figure  2.4  we  present  7  for  a  similar  general  body  in  polar  coordinates. 
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Figure  2.3:  Example  of  7:  general  body  in  Cartesian  coordinates. 

2. 3. 1.2  Extension  of  the  basis  function 

We  next  define  a  new  scalar  function  £!y  in  the  vicinity  of  T  which  represents  the 
numerical  counterpart  of  the  trace  of  the  solution  £p.  0,  is  obtained  by  extension 
of  from  T  to  7  using  Equation-Based  extension.  The  idea  behind  Equation- 
Based  extension  is  similar  to  the  approach  used  to  construct  compact  schemes  (see 
Section  1.2.1).  Given  the  function  and  the  first  normal  derivative  contained  in  £p, 
one  differentiates  the  Helmholtz  equation  Lq  u  =  F  to  obtain  higher  order  normal 
derivatives  and  uses  these  to  build  a  Taylor  formula  (2.24)  that  yields  the  value  of 
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£!y  at  every  node  7. 

v(n  +  5n,  s)  =  v(n,  s) +  (2.24) 

^  tl  orv 

We  denote  the  Equation-Based  Taylor  extension  as  ^  =  T(£g ,  F|r)  where  F 
is  the  source  term  of  the  equation.  An  example  of  Equation-Based  Taylor  extension 
in  polar  coordinates  can  be  found  in  Section  2.3.2  and  in  Section  2.3.5  we  show  the 
extension  in  elliptic  coordinates.  The  general  case  is  analyzed  in  [  ]. 

The  extension  T  applies  to  any  boundary  vector  function  =  {^Jq.  £,(/)-,  and 
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defines  a  new  scalar  function  in  the  vicinity  of  T  that  we  need  to  evaluate  at  the 
nodes  7.  The  operator  T(-,  ■,  0)  is  linear.  If  it  happens  that  ^  =  it|r  and  ^  =  |^|r, 
where  u  is  a  solution  to  the  Helmholtz  equation  Lq  u  =  F  on  Ilq.  then  T(£g,  Cq1 ,  T|r) 
approximates  u  at  the  nodes  7  with  the  accuracy  determined  by  the  order  of  the 
Taylor  formula  (2.24). 

The  number  of  terms  in  the  Taylor  formula  (2.24)  should  be  taken  as  the  min¬ 
imum  that  would  guarantee  the  design  rate  of  grid  convergence  of  the  overall  al¬ 
gorithm,  see  page  62  and  Section  Section  3.1.1.  Increasing  the  number  of  terms 
beyond  that  minimum  will  not  speed  up  the  convergence  any  further,  as  its  rate  is 
limited  by  the  order  of  accuracy  of  the  scheme.  Therefore,  in  practice  the  order  of 
the  Taylor  formula  (2.24)  is  always  kept  fixed.  In  theory,  however,  this  order  can  be 
allowed  to  increase.  Then,  by  invoking  a  Cauchy-Kowalevski  type  argument,  we  can 
conclude  that  the  series  will  converge  at  the  points  sufficiently  close  to  the  boundary 
as  long  as  all  the  data  are  analytic.  Among  other  classes  of  equations,  this  result 
holds  for  elliptic  PDEs,  which  is  our  case. 

Due  to  the  linearity  of  the  extension  T  one  extends  the  basis  functions  bn(s) 
in  expansion  (2.21).  For  an  interior  subproblem  we  define  1  (n)  =  T(bn,  0,0), 
£^1  (n)  =  T(0,  6n,0),  and  ^ ]  =  T(0,0,F|r),  where  the  latter  is  used  for  compu¬ 
tational  efficiency,  otherwise  the  contribution  of  the  source  term  is  unnecessarily 
calculated  several  times.  Thus  we  get 


£7,1  = 


(2.25a) 
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and 

Cjfi  =  cn7,1£ 7fi(n)>  (2.25b) 

and  we  define 

£7  =  £7,1  +  £74  +  £74-  (2.26) 

For  an  exterior  (homogeneous)  subproblem  we  define  £^j0(n)  =  T(bn,  0,0)  and 
£^o(n)  =  T(0,6n,0),  thus 

4  =  Ec«04(")  (2-27a) 

and 

tfffl  =  cn'°£"o(n)’  (2.27b) 

and  we  define 

$  =  £'o  +  £"0-  (2-28) 

Theoretically,  to  reach  nth  order  of  accuracy  for  the  overall  solution  to  the  prob¬ 
lem  one  requires  91  =  n  +  2,  see  [  ].  However  the  experimental  results  presented  in 

Section  3.1.1  and  in  [  ]  show  that  it  is  sufficient  to  take  91  =  n. 

We  should  also  mention  that  according  to  [  8,  27,  59],  for  maintaining  the  overall 
given  order  of  accuracy  across  the  domain,  it  may  be  sufficient  to  approximate  the 
boundary  conditions  with  lower  accuracy.  In  the  future,  it  may  be  of  interest  to 
investigate  whether  this  phenomenon  is  related  to  our  finding  that  the  order  of  the 
extension  operator  can  be  taken  lower  than  theoretically  predicted. 
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2. 3. 1.3  Difference  Potentials 

In  section  Section  2.1.3  we  required  an  arbitrary  function  w(x),  x  E  Mn,  that  satisfies 
the  Sonunerfeld  radiation  condition  (1.3)  and  such  that  it’s  vector  trace  is  equal  £p, 
i.e.  Trw  =  (w,wn)  =  (£/,£//)  =  £r  where  wn  is  the  normal  derivative  of  w  along 
r.  We  denote  the  discrete  trace  operator  similarly  as  the  continuous  one,  i.e.  Tr. 
In  the  discrete  case  the  trace  operator  become  the  value  of  the  function  at  the  grid 
boundary  7,  i.e.  Trw  =  w\^.  Thus,  we  require  the  discrete  w  to  satisfy  Trw  =  £7, 
where  £7  defined  in  Section  2. 3. 1.2. 

Numerically  one  creates  w(x)  using  the  same  structure  and  the  same  size  as  the 
solution  to  the  problem.  One  fills  w(x)  with  values  of  £7  at  nodes  corresponding  to 
7  and  zeroes  elsewhere.  We  present  such  a  procedure  in  Algorithm  1: 

Algorithm  1  Create  w  from  £7 

function  W(£7) 
w  E-  0  on  grid  N 

w\~i  ^  £7 

return  w 
end  function 

Let  <S(N,  k,  F,  g)  be  a  solver  of  the  equation  (1.1)  on  some  regular  grid  N  subject 
to  the  numerical  boundary  condition  w|grsj  =  g,  where  (9N  is  the  informal  notation 
for  boundary  nodes  of  N.  We  seek  to  solve  the  same  AP  with  different  source  terms. 
Therefore,  we  choose  a  fixed  g  to  provide  uniqueness  and  redefine  the  solver  as 
<S(N,  k,F).  Let  n  be  the  order  of  accuracy  of  S.  One  can  use  the  solver  described 
in  Section  1.2.1  or  as  described  in  more  detail  in  [  5,  12,  10,  55,  ].  We  denote  by 

S  the  solver  for  either  the  exterior  or  interior  problem. 


We  next  define,  in  Algorithm  2,  the  discrete  counterpart  of  Calderon’s  potential 
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Pnq  which,  unlike  in  continues  case,  will  take  as  input  w(x)  instead  of  £7  (see  Section 
2.1.3)  and  assume  that  Trw  =  £7  or  similarly  w(x) |7  =  £7,  i.e.  here  we  consider 
w(x )  is  known: 


Algorithm  2  Calderon  Potential  Pq9^7 

function  Pnq  (N,fc,'ir) 

>  Assuming  Trw  =  £7 

rhs  <-  0  on  grid  N 

rhs In,  <-  ( Lqw)\nq 

>  set  an  artificial  rhs 

v  <—  5(N,  k,  rhs ) 

return  w  —  v 

>  return  potential,  see  sections  2. 1,2. 3 

end  function 

( Lq  m)|Q9  is  the  result  of  the  direct  Helmholtz  operator  Lq  on  the  entire  grid  N  and 
truncated  to  the  domain  of  interest  Llq. 

One  defines  Calderon  Potential  in  the  following  way.  Let  u(s )  =  h(s)  and 
§^{s)  =  g(s)  therefore: 


Algorithm  3  Another  definition  of  Calderon  Potential 

function  Pnq(N,k,h,g,F) 

£7  T~{h,  3,  F ) 

w  <-  w(e7) 

return  Pq  (N,  k,  w) 

end  function 


The  difference  between  Algorithm  2  and  Algorithm  3  is  that  the  latter  creates  w  from 
u  and  using  an  extension  T  while  the  former  gets  w  as  an  argument.  Algorithm 
2  can  be  used  more  effectively  for  multiple  values  of  w,  therefore  Algorithm  3  is  not 
used  in  further  algorithms  and  presented  here  for  didactic  purpose  only. 


We  now  define  the  projection  operator,  which  numerically  become  a  truncation 


of  the  Calderon  operator  in  Algorithm  4. 
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Algorithm  4  Projection  operator 

function  P-?(N,fc,£7) 

w  <-  w(e7) 

return  Pnq(N,  k,  w)|7 

end  function 


2. 3. 1.4  Examples  of  solutions  to  different  problems 


We  next  solve  the  homogeneous  problem  Lqu  +  kqu  =  0  in  a  domain  Qq  (either 
interior  or  exterior)  and  let  the  boundary  condition  be  of  type  bcType  be  given  on 
dLlq.  We  assume  that  the  solution  has  an  expansion  in  some  functional  basis,  e.g. 
Fourier,  u{s)  =  and  f^O)  =  Endn(s)  and  §et 


Algorithm  5  The  Solver  for  a  homogeneous  problem 

function  Homogeneous-Solverq^  ( N,k,bcType,bc ) 

(£7,  Qi,  £7j,  Qii )  «-  BEP9(N,  k)  >  See  Algorithm  6 

(c7,  c11)  <—  Coefficients(Q/,  Qn,  bc-type,  be,  0)  >  See  Algorithm  7 

u  i —  0  on  grid  N 

n|n8^^(N,fc,^c7  +  ^V7))|n, 

return  u 
end  function 


Algorithm  6  used  in  Algorithm  5  plays  the  discrete  counterpart  of  (2.9)  which 
one  recasts  as  -P7£7  =  £7.  The  algorithm  Algorithm  6  doesn’t  provide  the  solution 
to  BEP,  but  the  matrix  with  columns  of  (P7  —  /)^7(n)  where  £7(n)  denotes  nth  basis 
function  extended  to  the  7,  see  Section  2. 3. 1.2.  This  matrix  is  used  to  construct  the 
linear  system  which  is  solved  to  find  the  coefficients  of  the  expansion  (2.21). 
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Algorithm  6  BEP 

function  BEPg(N,fc)  >  Calculate  operator  (P®  —  J)£7 

foreach  Basis  Functions  bn 

£7(n)  <—  T(bn ,  0, 0)  >  set  £7’s  as  columns  of  matrices  £7,  £7J 

£"(n)^T(0,&„,  0) 

5:  P/(n)  P®(N,  fc,£7(n))  >  solutions  to  AP’s  are  columns  of  P/,  P// 

Pjj(n)<-P*(N,fc,e?(n)) 

end  foreach 

Qi  <-  Pi-  £7 

Q//  t—  P/j  —  ^ 

10:  return  (£*,  Qj,  Qjj) 

end  function 


The  coefficients  are  actually  found  using  Algorithm  7  in  the  least  square  sense  using 
QR,  since  the  linear  system  is  redundant  for  sufficiently  fine  grid.  It  is  solved  by 
least  squares  via  QR.  See  [42]  and  following  examples. 


Algorithm  7  Compute  coefficients  of  an  expansion  on  an  interface 

function  Coefficients (Qi,Qn,bc-type,bc,InHomoPart) 

c  <—  CoeffsOf(foc)  >  coefficients  of  be  =  J2n  cnbn 

switch  be  Type 

case  Dirichlet  >  «|7  =  be 

5:  c1  c 

c11  «—  QR(Q//,  —Qi  Ci  —  InHomoPart ) 
case  Neumann  >  rtn|7  =  be 

C11  <r~  C 

cl  <r-  QR(Q/,  —Qn  ci i  —  InHomoPart) 

10:  case  Robin  >  (era  +  /3un)|7  =  be 

c1  <—  QR  (Qi  —  |  Qn ,  —  Qn  c  —  InHomoPart) 

C11  7-  i 0-  “cj 
c  ^  /3C 

end  switch 

return  (cz,  c11) 

15:  end  function 


We  next  solve  the  inhomogeneous  problem  in  Algorithm  8.  A  distinctive  feature  of 
compact  schemes  is  that  the  right-hand  side  of  the  difference  equation  gets  trans¬ 
formed,  see  (1.8),  we  denote  this  transformation  as  B(F).  However  the  transforma- 
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tion  of  the  source  term  requires  F  to  be  defined  not  only  on  Qq  where  it  is  originally 
defined  in  the  analytical  problem,  but  also  on  7.  Therefore,  we  continue  F  from  Qq 
to  Qq  U  7  using  the  Taylor  extension  T(F). 


Algorithm  8  The  Solver  for  an  inhomogeneous  problem 

function  Inhomogeneous-Solver^  (N ,k ,bcType,bc) 

BEP9(N,fc) 

<—  Ex( 0, 0,  F)  >  calculate  inhomogeneous  part  of  Ex  only  once 

Pf  4—  P^{ N,  >  which  give  another  AP  to  solve 

5:  Qf  4-  Pf  — 

GF  =  S{N,k,B(T(F))) 

(c1 ,  c11)  4—  Coefficients^/,  Qn,  be  Type,  be,  Qf  +  GF\ 7) 
u  4 —  0  on  grid  N 

10:  return  u  +  GF 

end  function 


We  finally  solve  the  transmission-reflection  problem  driven  by  the  incident  wave 
u(mc)  jn  Algorithm  9.  The  media  of  the  interior  part  can  be  heterogeneous.  The 
exterior  problem  is  considered  in  the  far  field  and  therefore  the  media  is  homoge¬ 
neous.  However,  a  jump  is  allowed  at  the  interface  between  the  interior  and  exterior 
problem.  The  solution  to  the  exterior  problem  is  considered  as  a  sum  of  the  scattered 
and  incident  field. 

A  clear  advantage  of  our  method  is  that  scattering  about  a  given  shape  but 
for  multiple  angles  of  incidence,  and  even  for  different  boundary  conditions,  can 
be  computed  very  efficiently.  This  is  particularly  important  if  the  direct  scattering 
problem  needs  to  be  solved  many  times  while  using  an  iterative  method  to  solve  an 
inverse  scattering  problem. 
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Algorithm  9  The  Solver  for  Transmission-Reflection  problem 

function  TRANSMISSION-REFLECTION-SOLVER(N,fc0,fcl^cJ?/pe,Z^lnc^) 

(#J,  Qi,i,  Qi,n)  <-  BEPiCNi,  h) 

(£’I,Qo,i,£’II,Qo,n)  BEPo(No,fco) 

£  r  ^ 


12: 


£n  r 


( ^ 

6:  £!•*’<- £7a:(0,  0, -F1) 


Qo,(inc)  t— 

GFi  =51(N1)&i,B(T(F))) 
Qo,i,  Qoji 
Qi,i,  Qi,ii 

u  ■< —  0  on  entire  grid  No  U  Ni 

u\n0  ^  S°(Ho,ko,£c)\n0 


-  U^nc) 


c  =  QR 


Qo,(inc) 

— GFi\y  —  Qi,e 


M(inC)|on 


w|f2i  <—  S\l 

return  u 
end  function 


,ki,£c  +  £F)\n1+GF1 


>  Interior 
>  Exterior 


Consider  a  set  of  incident  waves  u^mc\6m),  The  naive  algorithm  uses  Algorithm 
9  m  times. 


Algorithm  10  Inefficient  Solver  for  Transmission- Reflection  problem  with  multiple 
incident  angles 

function  Transmission- REFLECTiON-SoLVER-MA(N,fco^i,&cJ2/pe,'u(inc\0) 
foreach  Incident  Angle  9m 

Um  <—  Transmission  —  Reflection  —  Solver{  N,  k,  be -type,  ut™^) 

end  foreach 
return  U 
6:  end  function 


However,  it  can  be  done  much  more  efficiently,  since  the  call  to  the  most  time 
expensive  function  BEP  does  not  depend  on  the  impinging  wave  and  therefore  can 
be  called  once  per  problem  (Interior,  Exterior)  regardless  of  the  number  of  incident 
angles  used: 
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Algorithm  11  Efficient  Solver  for  Transmission- Reflection  problem  with  multiple 
incident  angles 


function  TRANSMISSION-REFLECTION-SOLVER-MA(N,fco,/ci,&C_h/pe,u(inc),0) 

{Q1  iQuiQ11  ,Qi,ii)  BEPi(Ni,fci)  >  Interior 

(£°’J,  Qo,i,  Zy’11,  Qo,ii)  BEP0(N0,  ko)  >  Exterior 

£7  «-  (  £7’J,  Q11  V 
$  <-  (  #T,  f 

6:  QF  <-  Ex(0, 0,  F) 

Qi,F<-P$(tii,ku£l'F)-$F 

GF1  =S1(N1,k1,B(T(F))) 


foreach  Incident  Angle  0m 

Q0,(inc)  P7°(N0,fc0,^nC)|7)  -  ^C)|7 


c  =  QR 


Qo,i,  Qoji 
Qij,  Qiji 

12:  u  <—  0  on  entire  grid  No  U  Ni 

R|n0  <5°  (No,  ko,  C7c)|o0  +  u<rn  ^|fl0 
u\ni  ^S1(N1,k1,Qc  +  QF)\Qi+GF1 


Q  0,(inc) 

—GF\\y  —  Qi,f 


Um  U 

end  foreach 
return  U 
18:  end  function 


2.3.2  Cartesian  Coordinates:  Homogeneous  Dirichlet  Problem  in 
a  Circle 

Consider  the  Dirichlet  problem  (2.29)  in  a  circular  domain  Pi  =  {(x,  y)\ \J  x2  +  y2  < 
R}  solved  on  a  Cartesian  grid,  where  the  boundary  is  not  aligned  to  the  grid. 

A u  +  k2u  =  0  x2  +  y2  <  R2, 

<  (2.29) 

u(x,  y)  =  g(x,  y)  x2  +  y2  =  R2. 

We  introduce  a  finite-difference  approximation  to  this  problem.  Let  N  be  the 
Cartesian  grid  on  the  (—xt,  xt)  x  (—yt,  yt),  where  xt  =  yt  >  R  and  let  pij  =  ( Xi,yj )  E 
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N  be  the  nodes  of  the  grid  with  \pij\  =  w  x'f  +  y1- .  Denote  M  =  {pij  :  \pij\  <  R}. 


Figure  2.5:  Circular  body  in  Cartesian  coordinates  with  7  defined  in  (2.23). 

To  solve  problem  (2.29)  with  higher  order  accuracy  we  first  define  the  grid  coun¬ 
terpart  of  original  boundary  shape  T  whose  fragment  and  approximation  are  shown 
in  Figure  2.5  as 

7  =  N+nr,  (2.30) 

where  is  the  9-point  stencil  centered  at  node  pl.j.  see  Figure  2.2.  N+,N“  are 
non-empty  sets  defined  as 

N+=  U  np«.  U  fw 

Pi.jSM  Pi, jSN\M 
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2.3.2. 1  Equation  Based  Extension 

The  set  7  introduced  by  formula  (2.23)  will  be  used  to  define  the  density  £7  of  a 
difference  potential.  In  turn,  £7  will  be  obtained  from  by  means  of  the  Equation- 
Based  Taylor  extension  as  explained  in  Section  2. 3. 1.2.  We  present  an  example 
of  a  5th  order  Equation-Based  Taylor  extension  (2.31).  Since  the  boundary  shape 
is  a  circle  it  is  convenient  to  rewrite  the  boundary  condition  as  u{R,9)  =  g{9). 
Although  we  solve  the  homogeneous  problem,  the  extension  is  for  a  more  general 
inhomogeneous  case  of  an  equation  based  extension. 

Sr3  Sr 4 

u(R  +  Sr,  9)  =  u(R,  9)  +  Srur  H — —urr  -| — —urrr  +  urrrr .  (2-31) 

We  use  Helmholtz  equation  in  polar  coordinates 

urr  +  -Ur  +  ~2u0e  +  k2{r,  9)u  =  F(r,  9)  (2.32) 

and  solve  (2.32)  for  the  second  radial  derivative 

P  1  1  *2 

Urr  =  r - Ur - ^Uqq  —  k  U. 


We  differentiate  it  twice,  for  the  third  derivative 
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and  the  fourth  derivative 


ur 


„  (  2  1  2\/6  4  1 

—  L rr  4“  I  —Urrr  —^Ur  I  A  I  —jUgg  4  ^Urgg  ^)Urrgg 

\r  r  r°  J  \  r  r°  r 

—  k2urr  —  Akkrur  —  2  k2u  —  2  kkrru 


=  Frr  —  2  (k2  +  kkrr)  u  —  2  ^ ^  +  2 kkr^j  ur  +  —  k2^j 


ur 


6 

r4 


4 

—  Urrr  —r  Uqq  ~\  o  Ur 00  “ n'U'OOrv 

r*^ 


The  mixed  derivative  Uggrr  is  obtained  using  the  angular  derivative  of  the  equa¬ 
tion: 

w  1  1  ,2 

Uggrr  =  rgg  —  ~Urgg  —  —^Ugggg  —  K  Ugg. 

Finally, we  define 


7^ Tru  =  T{u,ur,F,\pij\-  R),  pid  <E  7, 


where  T  denotes  the  Taylor  extension  (2.31). 


2. 3. 2. 2  Auxiliary  Problem 


Consider  the  following  auxiliary  problem 


LpijUpij  —  fpi,j  Pi,j  £  N  \ 
l. Pi,jVPi,j  =  0  Pi,j  £  <9F1> 


(2.33) 
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where  the  boundary  condition  lpt  jvPi  .  =  0  in  (2.33)  can  be  arbitrary  as  long  as  it 
guarantees  the  existence  and  uniqueness  of  the  solution  for  any  right-hand  side  fPi  . , 
and  well-posedness  of  the  entire  formulation.  In  practice,  the  AP  is  also  chosen  so 
that  it  admits  an  easy  numerical  solution.  For  example,  one  can  use  an  absorbing- 
type  boundary  condition  in  (2.33): 


L  v  =  < 


(vx)nj  =  ikvPij 


vPi,j  —  o 


xi,j  e  {-xt,xt}, 

y%,i  e  {-vuvt}- 


(2.34) 


The  AP  (2.33)  is  used  to  compute  the  difference  potential  with  the  density 
defined  on  the  grid  boundary  7.  For  a  given  £7,  we  first  introduce  the  auxiliary 
function 


£7 1  Pij  PiJ  ^  7) 

0  elsewhere, 


(2.35) 


so  that  TrwPi  .  =  ^7 .  Then  we  define  the  right-hand  side  for  the  AP  as  follows 


fp%,j  ~  " 


Lpi,jwPij  Pi,j  ^ 

0  elsewhere. 


(2.36) 


2. 3. 2. 3  Reduction  To  The  Boundary 

Denote  by  vPi  .  the  solution  to  the  auxiliary  problem  (2.33)  with  the  right-hand  side 
given  by  (2.36).  The  difference  potential  with  the  density  £7  is  defined  for  the  grid 
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nodes  pij  E  N+  according  to  the  following  formula: 

=  WPi,j  ~  VPi,ji  (2.37) 

where  wPi  -  is  given  by  (2.35).  The  difference  potential  (2.37)  converges  to  the 
continuous  potential  (2.16)  with  the  design  rate  as  the  grid  size  decreases,  see  [  , 

and  also  [  ]. 

We  next  assume  the  Fourier  expansions,  i.e.  bn(0)  =  el™ ,  see  Section  2.3. 1.2: 

M 

u(R,  9)  =  (2.38a) 

n=—M 

and 

M 

ur(R,6)  =  ^2  ur[n\em8 .  (2.38b) 

n=—M 

When  we  solve  the  Dirichlet  problem  the  coefficients  u[n]  are  known  while  ur[n] 
are  unknowns.  We  next  define  two  sets  of  the  source  term  via  £7: 

=  %j(ein0,O)  =  T(ein6,  0,  0,  <5r),  n  E  [— M,  M\ 

and  another  for 

e57( n )  =  77j(0,  eine)  =  T(0,ein9,0,6r),  n  E  [— M,  M\, 

where  <5r  represents  the  shortest  distance  between  pij  E  7  and  the  analytical  shape 
r,  i.e.  in  this  instance  hr  =  \Pi,j\  ~  R- 
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We  next  introduce  the  discrete  projection 


Py£,y  — 


and  the  difference  counterpart  of  BEP  (2.9): 


—  £7 


or  equivalently 


(wPi,j  vPij)  I7  —  wPi,j\y 


We  rewrite  the  BEP  (2.39)  as 


Py£y  £7  —  0> 


and  define 


Qq€  7  —  (■^>?  -0£y  —  (wPij  vPi.j')\i  wPi,j  It  —  vPi,j\y 


We  next  use  (2.26)  to  recast  (2.40)  as 


Qq€ 7  —  Qq£y  +  Qq^y  , 


(2.39) 


(2.40) 


where 
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see  (2.25).  We  next  combine  Qq£y(n)  and  Qql,-/ (n)  in  a  matrix: 


/  | 

1 

1 

\ 

Q  =  [QAQii]  = 

Qq$(-M)  • 

•  Qq^(M) 

Qgtfi-M)  • 

•  Qqe/(M) 

V  1 

1 

1 

1  / 

(2.41) 

The  size  of  the  matrix  Q  is  I'y |  x  2(2 M  +  1),  where  |,y|  denotes  number  of  nodes  of 
grid  in  7,  which  also  the  size  of  Qq£,™(n)  and  £™(n).  Thus,  we  are  looking  for  a 
non-zero  vector 


c  =  [u[-M], . . .  ,u[M],ur[-M\, . . .  ,ur[M]]T  (2.42) 

that  satisfies 

Qc  =  0.  (2.43) 


2. 3. 2. 4  The  Solution 

Equation  (2.43)  has  multiple  solutions  since  we  did  not  yet  use  the  boundary  con¬ 
dition,  which  we  do  now.  The  problem  (2.29)  is  a  Dirichlet  problem  and  therefore 
the  2 M  +  1  coefficients  of  (2.38a)  are  easy  to  calculate  by  Fourier  transforming  the 
data  g(0).  One  then  denotes 

c1  =  [u[— M], . . . ,  u[M]]T , 

c11  =  [ur[-M},...,ur[M}}T. 


(2.44) 

(2.45) 
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and  rewrite  (2.43)  as  a  system  of  equations  with  respect  to  c11 

Qnc11  =  -QIcI.  (2.46) 

For  a  sufficiently  fine  grid,  system  (2.46)  is  over  determined,  and  its  solution 
is  to  be  sought  for  in  the  sense  of  the  least  squares.  However,  as  the  solution  to 
the  original  continuous  problem  exists,  and  system  (2.46)  equivalently  represents  its 
fourth  order  accurate  finite  difference  approximation,  the  residual  of  its  least  squares 
solution  is  expected  to  vanish  with  the  rate  0(/i4)  as  the  grid  is  refined.  In  this  sense, 
the  overdetermined  system  (2.46)  can  be  said  to  have  an  “almost  classical”  solution. 

To  obtain  the  solution  to  problem  (2.29)  one  again  solves  the  auxiliary  problem 
(2.33).  Thus,  the  approximation  is  given  by  u  =  Pq^ 7  where 

£7  = 

n  n 

where  now  all  the  coefficients  are  known. 

2.3.3  Cartesian  Coordinates:  Homogeneous  Neumann  Problem  in 
a  Circle 

Consider  the  Neumann  problem: 

Au  +  k2u  =  0  x2  +  y2  <  R2, 

< 


ur(x,  y)  =  g(x,y)  x2  +  y2  =  R2. 


(2.47) 
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The  solution  to  the  problem  (2.47)  is  similar  to  the  solution  of  (2.29).  The  main 
difference  is  that  unlike  in  (2.29)  where  first  2 M  +  1  elements  of  c  are  known  in 
Neumann  problem  the  second  2 M  +  1  elements  of  c  are  known  and  the  first  2 M  +  1 
elements  are  unknown.  Thus,  to  solve  the  Neumann  problem  one  exchanges  between 
the  right  hand  side  and  left  hand  side  of  (2.46),  i.e. 

Qic1  =  -Q,,cn  (2.48) 

and  the  rest  of  the  procedure  remains  unchanged.  Therefore,  one  computes  un¬ 
known  coefficients  using  (2.48)  and  approximates  the  solution  to  problem  (2.47)  by 
computing 

=  Y1  “NC5(n)  +  («) 

n  n 

and  then  solving  an  auxiliary  problem  (2.33)  for  u  =  Pq  H ■ 

2.3.4  Cartesian  Coordinates:  Homogeneous  Robin  Problem  in  a 
Circle 

We  next  consider  the  Robin  problem: 

A u  +  k2u  =  0  x2  +  y2  <  R2, 

<  (2.49) 

om(x,  y )  +  (3ur(x,  y)  =  g(x,  y)  x2  +  y2  =  R2, 

where  a2  +  j32  7^  0.  One  uses  the  orthogonality  of  the  basis  in  the  expansions  (2.38) 
to  rewrite  the  boundary  condition  as  au[n ]  +  /3ur[n]  =  g[n],  Vn  G  [— M,  M].  Hence, 
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one  solves  the  following  system  of  equations 

ac1  +  /3c11  =  g,  (2.50) 

Qic1  +  Qnc11  =  0.  (2.51) 


One  solves  (2.50)  to  find  c 


ii 


ii  1  -  a  i 

c  =~P~r 


and  then  uses  it  in  (2.51)  to  obtain 


Qic1  +  Qid^g-  j^c1)  =  0 


and  hence  one  solves 

[  Qj  -  jj  Qii)  c1  =  ~  -jj  Qii9 

to  find  c1  and  substitute  it  in  (2.52)  to  obtain  c11 . 
Alternatively,  one  solves 


( Qii  -  —Qi)cn  =  — Qig 
a  a 

to  compute  c11 ,  and  then  substitute  it  in  c1  =  -g  —  ^c11 . 
Once  the  coefficients  c1 ,  c11  are  known,  one  then  computes 


£7  =  ^2u[n]Cif(n)  +  ^2ur[n}^I(n) 

n  n 


(2.52) 


(2.53) 


to  approximate  the  solution  to  the  problem  (2.49)  by  solving  an  auxiliary  problem 
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(2.33)  for  u  =  Pq£7. 


2.3.5  Cartesian  Coordinates:  Inhomogeneous  Equation  in  an  El¬ 
lipse 


We  next  present  the  inhomogeneous  Dirichlet  problem  (2.54)  for  an  elliptical  body 
to  be  solved  on  a  Cartesian  grid,  see  Figure  2.6. 

An  +  k2u  =  F(x,  y)  p-  +  fx  <  1> 

(2.54) 

2  2 

u(x,y)  =  g(x,y)  p  +  p  =  1, 

where  a  =  d cosh iio,b  =  dsinhr/o  are  the  major  and  minor  semiaxes  of  the  ellipse 
respectively,  d  =  V a2  —  b2  is  the  focal  distance  and  r/o  is  the  elliptical  radial  coordi¬ 
nate.  One  reexpresses  (2.54)  using 


??(x,  y)  =  Re 


/  x  +  iy\ 

^arcosh  -  ) 


as 

An  +  k2u  =  F(x,  y)  p(x,  y)  <  rjo, 

<  '  (2.55) 

u{x,y)  =  g(x,y)  rj(x,y)  =  rj0- 

The  finite-difference  approximation  to  this  problem  is  given  by  an  auxiliary  prob¬ 
lem  (2.33).  The  grid  N  is  the  Cartesian  grid  on  (—xt,xt)  x  (~yt,Ut),  see  Figure 
2.6  and  pij  =  ( Xi,yj )  G  N  is  the  node  of  the  grid.  Instead  of  \pi,j\  we  define 
rjij  =  p{pij)  =  rj(xi,yj)  and  denote  M  =  {pij  :  r/jj  <  rjo}-  The  definition  of  7 
remains  unchanged. 
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Figure  2.6:  Example  of  discrete  auxiliary  problem:  elliptical  body  in  Cartesian  coordinates. 


Obviously,  since  the  boundary  is  not  aligned  to  the  grid  one  needs  to  develop  an 
Equation  Based  Taylor  extension  for  the  new  curve.  The  general  case  is  explained  in 
] .  Here  we  limit  our  discussion  to  the  Helmholtz  equation  in  elliptical  coordinates 


h2{r],p) 


(uvv  +  uw)  +  fe2(?y,  p)u  =  F(rj,  p), 


(2.56) 


where  h(rf,  p)  =  dy/ sinh2  7]  +  sin2  p  is  the  scale  factor. 


We  first  solve  (2.56)  for  uvv  to  get 


uvv  =  h2F  —  uw  —  h2k2u 


then  we  differentiate  it  to  get  the  third  derivative 

uvrir)  =  2 hhvF  +  h2Fv  —  —  2  ( hhvk 2  +  h2kkv )  u  —  h2k2uv 


85 


82  CHAPTER  2.  PROBLEMS  WITH  NON-ALIGNED  INTERFACES 

and  the  fourth  derivative 


v^rjrjr]7]  —  2  (/r“  +  hh^)  F  +  ihhr] Eri  +  h 

—  2  ((h2  +  hhvv)  k 2  +  Ahkh^kr,  +  h2fcfcw  +  h2k 2)  u 

—  4  ( hhvk 2  +  h2kkv)  uv  —  h2k2um 


and  similarly 

—  2  (/i^,  T  hh<p<p^  F  T  IhhpFp  -)-  h  Ecp<p 

—  2  ( (/i2  +  hhpp)  k2  +  Ihkhpkp  +  h2kkvip  +  h2k 2 )  u 

—  4  ( hhcpk 2  +  h2kkv)  u ^  —  h2k2UtpV 


We  next  use  these  derivatives  in 


u(r}  +  h,  ip) 


h 2  h3  h 4 

4~  kurj  T  Ur/ri  H  T  irnvn 


(2.57) 


and  define 

7^°  Tru  =  T{u,  uv,  F,  rji:j  -rj0),  pij  E  7, 


(2.58) 


where  T  denotes  Taylor  extension  (2.57). 


We  next  assume  the  Fourier  expansions 

M 

u(r]0,ip)=  ^  u[n\eirvp 

n——M 


(2.59a) 
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uv(vo,T)=  ^rj[n]emv>  (2.59b) 

n=—M 

and  define  the  sets  of  £™(n).  In  order  to  reduce  the  computational  cost  instead  of 
two  sets 

e{(n)  =  T(einx,  0,  F,  St),  n  G  [-M,  M } 

and 

^(n)  =  T(0,  einx,  F,  Sr),  n  G  [-M,  M], 

where  <5r  =  ViJ  ~  Vo  is  the  shortest  distance  from  pij  6  7  to  T.  One  solves  these 
two  sets  without  the  source  term  i.e. 

^(n)=T(einx,0,0,6r),  n€[-M,M } 

and 

^(n)  =  T(0,  einx,  0,  Sr),  n  G  [-M,  M } 

and  defines  an  additional  AP  for  the  source  term  F.  This  is  done  so  we  need  to 
compute  it’s  part  of  the  extension  only  once 

$  =  T(0, 0,  F,  Sr),  n  G  [-M,  M\. 

This  adds  another  column  to  the  matrix  Q  in  equation  (2.43),  Qq^  =  (Pq  -  I)tf. 
More  precisely  there  is  another  change  to  (2.43),  i.e.  it’s  become  inhomogeneous: 


Qc  =  —GF  |7, 


(2.60) 
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where  GF  is  the  solution  of  AP  (2.33)  with 


•f o ! .  j 


BrFPiJ 


Pij  e  M  U  7, 
elsewhere , 


(2.61) 


where  the  operator  B  represents  the  right  hand  side  stencil  of  the  compact  scheme, 
see  Section  1.2.1,  particularly  the  right  hand  side  of  (1.8).  Note  that  B  requires  F 
to  be  defined  for  pij  for  which  r]ij  >  ?yo  where  F  is  not  necessary  known.  Therefore, 
one  extends  the  source  term  F  to  these  nodes  using  regular  (non  equation  based) 
Taylor  extension  which  we  denote  T.  See  also  [  ]. 

Thus,  the  equation  (2.46)  becomes 


Qiiu^[n]  =  -  Qiu[n]  -  Qqtf  -  GF |7.  (2.62) 


To  obtain  the  solution  to  problem  (2.55)  one  computes  the  coefficients  from 
(2.62)  and  then  solves  the  auxiliary  problem  Pq  using 

£7  =  +  ^^[n]^1  (n)  +  Qgtf. 

n  n 


Finally  u  =  Pq  £7  +  GF. 
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2.3.6  Polar  Coordinates:  Scattering  about  an  Ellipse 


Consider  the  scattering  problem  in  polar  coordinates  about  an  elliptical  body 


usrr  +  \usr  +  usee  +  k2us  =  0  rj(R,  8)  >  rj0 


us(R,  9)  =  g(6) 


r)(R,9)  =  77o, 


(2.63) 


w  r— oo 


where  us  denotes  a  scattering  field,  g{0)  =  —  u^mc\d)  =  e  lk(xco sO+ysme)  jg  an 
incident  wave,  7/q  is  the  elliptical  radial  coordinate  and 


and  d  is  the  semi-focal  distance,  see  Figure  1.2. 


Since  on  a  computer  one  can’t  solve  the  problem  in  an  infinite  domain  one  trun¬ 


cates  the  infinite  domain  and  defines  an  artificial  boundary.  It  is  convenient  to  choose 
the  artificial  boundary  to  be  aligned  to  the  the  grid,  therefore,  the  numerical  domain 


becomes  a  ring,  Figure  2.7.  The  Sommerfeld  condition  lim  r  2  (usr  +  ikus)  =  0  is 


r— >00 


changed  to  an  Absorbing  Boundary  Condition  (ABC)  which  uses  the  decomposition 
of  the  waves  into  incoming  and  outgoing  and  allows  the  propagation  of  only  the 
outgoing  waves.  Such  an  approach  resembles  the  idea  of  Sommerfeld  condition  of 
“no  wave  radiating  from  infinity”  by  “no  wave  is  entering  the  computational  do¬ 
main”  or  more  precisely  no  wave  is  reflecting  from  the  artihcial  boundary.  This  is  a 
commonly  used  approach  [  1,  29],  and  is  valid  since  we  are  considering  only  regular 
bounded  solutions,  which  eliminates  the  possible  ambiguity  in  the  definition  of  the 


89 


86  CHAPTER  2.  PROBLEMS  WITH  NON-ALIGNED  INTERFACES 


Figure  2.7:  Example  of  discrete  auxiliary  problems:  elliptical  body  in  polar  coordinates. 

incoming  and  outgoing  waves  pointed  out  in  [  ]. 

Occasionally,  for  an  elliptical  or  oval  like  scatterer  one  wishes  that  the  outer 
artificial  surface  resemble  the  scatterer  to  prevent  unnecessary  interior  nodes.  One 
then  uses  an  ABC  in  elliptical  coordinates  (see  [  ,  43,  1  ]).  However,  we  set  the 

artificial  boundary  on  an  outer  circle  of  the  grid.  We  Fourier  transform  the  auxiliary 
problem  (2.33)  to  get  a  linearly  solved  tridiagonal  linear  system  [  ],  see  algorithm 

described  in  [  ].  Therefore,  in  this  instance  we  neglect  the  problem  of  redundant 
nodes. 


The  AP  (2.33),  transformed  to  the  Fourier  space,  has  the  advantage  that  the 
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ABC  becomes  an  exact  condition  instead  of  an  approximation  of  the  Sommerfeld 
condition  [  ] .  The  boundary  condition  on  the  inner  circle  may  be  chosen  arbitrarily 
since  unlike  the  ABC  it  belongs  to  the  AP  but  not  to  the  original  problem.  One  de¬ 
fines  a  homogeneous  boundary  condition  on  the  inner  circle.  By  this  we  accomplish 
the  new  definition  of  lPij  in  (2.33). 

One  then  uses  an  Equation  Based  Taylor  extension  7^°  defined  in  (2.58)  (see 
]  for  more  general  case)  and  assumes  the  Fourier  expansion  of  the  solution  on  an 
interface  T  as  given  in  equations  (2.59a)  and  (2.59b).  Next,  one  solves  the  AP  for 

£y(n)  =  T(einx,0,0,dr),  ne[-M,M } 


and 

e^(n)  =  T(0,  einx,  0,  <5r),  n  G  [— M,  M], 
where  hp  =  Tji.j  —  r/o  is  the  shortest  distance  from  pij  S  7  to  T. 

Finally,  one  defines  the  matrix  Q  and  using  (2.40),  one  solves 

Qiiu^n]  =  (2-64) 

where  us  =  — n(mc)  and  computes 

£7  =  ^2us[n]^(n)  +  J^[n]^7(n) 

n  n 

to  approximate  the  solution  to  problem  (2.63)  as  us  =  Pq  £7  or  the  total  field  as 
u  =  us  +  u (inc)  =  Pq  ^  +  u(inc) . 
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2.3.7  Transmission— Reflection  problem 


Consider  the  problem  (1.2)  described  also  in  Figure  1.1.  For  convenience  we  redefine 
it  as 

A u  +  k%u  =  0  a:  G  Oo, 

<  (2.65) 

Au  +  k\(x)2u  =  F(x)  x  6  fii. 


AUX-BC 


Figure  2.8:  Interior  and  Exterior  subproblems 
In  this  instance  we  consider  that  T  is  an  ellipse  and  therefore  the  auxiliary 
problems  that  solve  system  of  BEPs  (2.17)  becomes  the  problems  (2.54)  and  (2.63). 
When  these  problems  are  solved  the  information  about  the  boundary  function  of 
these  problems,  g,  are  used  at  a  very  late  stage.  More  precisely  the  matrices  Q  in 
both  problems  are  computed  without  any  knowledge  about  g  and  this  fact  is  a  key 
property. 

Thus,  one  defines  four  sets  of  £7,  two  sets  (n),  £7fo(n)  f°r  exterior  and  two 
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sets  ^  i(n),  £yfi(n)  for  the  interior  subproblem,  i.e.  assume  an  expansion  (2.59a), 
(2.59b)  and  obtain  a  matrix  Q. 

Let  Q 1  be  the  matrix  Q  of  the  interior  subproblem  (2.54)  for  which  the  BEP  is 
given  by 

Q1c  =  -Q1$- 

and  let  Q°  be  the  matrix  Q  of  the  interior  subproblem  (2.63)  with  it’s  BEP 

Q°c  =  Qo4inc)- 


When  problems  (2.54)  and  (2.63)  are  solved  half  of  the  vector  of  coefficients  c  are 
known  from  the  boundary  condition  and  the  other  half  are  solved.  In  this  instance 
c  is  fully  unknown,  and  instead  one  solves 


'  Q'  )  =  (  Vi  ft  -  Gl„ |7  \ 

v  0"  j  C  ”  {  Qof?nc|  j 


(2.66) 


which  is  again  overdetermined  for  a  fine  enough  grid  and  it’s  solution  exists  again 
as  long  as  solution  to  original  problem  does.  Hence  (2.66)  is  understood  in  the  least 
square  sense. 


To  approximate  the  solution  to  problem  (2.65)  one  computes 


4  =  ^u[n]£!y{n)+^u^[n]£I1I(n)  +  Q1tf. 


n 


n 


93 


90  CHAPTER  2.  PROBLEMS  WITH  NON-ALIGNED  INTERFACES 


and 


and  solves 


£7  =  +  2SM£57(«) 

n  n 


Px  e( + gxf 

Po  [£°  -  4inc)]  +  «(inc) 


x  S  Qi, 


x  £  Oq  - 


2.4  Complexity 

The  key  contribution  to  the  overall  complexity  of  the  proposed  algorithm  is  the 
repeated  solution  of  the  discrete  AP  (2.33).  Each  subproblem  of  the  transmission 
reflection  problem  needs  to  be  solved  (2 M  +  1)  times  to  find  the  coefficients  and  one 
more  to  obtain  the  solution  when  coefficients  are  known.  Altogether,  it  needs  to  be 
solved  4 (M  +  1)  times.  However,  only  the  right-hand  side  of  the  AP  changes  from 
time  to  time,  whereas  the  two  operators  remain  the  same. 

In  the  case  of  constant  coefficients,  each  of  the  4(M+1)  equations  requires  a  FFT 
solve  which  has  a  log-linear  complexity  with  respect  to  |N|,  the  dimension  of  the 
problem,  i.e.  the  number  of  nodes  of  the  grid  N.  In  the  case  of  variable  coefficients, 
the  overall  complexity  will  be  that  of  a  single  sparse  LU  decomposition  of  a  matrix 
of  dimension  |N|  x  |N|  plus  4 (M  +  1)  backward  substitutions.  This  is,  of  course, 
significantly  faster  than  solving  the  overall  system  many  times. 

In  the  case  of  variable  coefficients,  one  can  use  an  iterative  solver,  for  example 
see  [  ].  In  3D,  a  straightforward  Gaussian  elimination  (LU  decomposition)  is  not 

feasible.  For  constant  coefficients,  a  Fourier  based  solver  that  does  not  require 
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storing  the  matrix  and  has  a  log-linear  complexity  will  provide  the  most  efficient 
approach  to  solution,  if,  of  course,  the  AP  can  be  formulated  so  that  it  will  admit 
the  solution  by  FFT.  Otherwise,  an  iterative  solve  will  become  a  necessity  in  3D 
for  either  variable  or  constant  coefficients.  For  example,  using  the  iterative  scheme 
Risolv  of  [  ],  one  needs  to  find  the  optimal  parameters  of  the  algorithm  only  once 

independent  of  the  RHS.  Let  Nnz  be  the  number  of  non-zero  entries  in  the  system 
matrix  that  has  the  dimension  |N|  x  |N| ,  Nk  be  the  dimension  of  the  Krylov  subspace, 
and  A'a  be  the  number  of  times  we  apply  the  Arnoldi  algorithm.  Then  the  total 
amount  of  work  for  the  first  solve,  i.e. ,  for  one  RHS,  is  approximately  A^x  A/kX  Nnz+ 
N\  x  \N^  x  |N| .  For  each  subsequent  solve  the  total  work  is  approximately  N\  x  Nk  x 
Nnz.  Hence,  it  no  longer  depends  on  A^|N|,  which  means  that  the  dependence  on  the 
dimension  |N|  of  the  system  matrix  disappears,  and  the  dependence  on  the  dimension 
Nk  of  the  Krylov  subspace  becomes  linear  rather  than  quadratic.  Regarding  the 
preconditioning  of  the  Helmholtz  equation,  see,  e.g.  [  1,  19,  IS,  20]. 

Another  contribution  to  the  overall  complexity  is  the  QR  decomposition.  If  the 
modified  Gram-Schmidt  algorithm  is  used,  then  the  corresponding  cost  is  about 
2(2 M  + 1)2|7|  operations,  where  for  the  current  two-dimensional  setting  |q|  f-KJ 
The  cost  of  all  other  components  of  the  algorithm,  see  Section  2.3.1,  is  negligible. 

We  re-emphasize  that  as  the  discretization  grid  is  refined,  only  the  quantities 
|N|  and  |y|  increase,  whereas  M  stays  the  same.  This  is  the  strategy  we  have 
adopted  for  all  our  numerical  simulations,  see  Chapter  3.  Of  course,  choosing  M 
grid-independent  represents  the  most  conservative  scenario.  It  allows  us  to  make 
sure  that  the  accuracy  of  the  spectral  representation  at  the  boundary  (in  the  basis 
of  dimension  M)  will  certainly  exceed  any  accuracy  that  one  may  possibly  obtain 
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on  the  grid.  In  fact,  M  does  not  always  have  to  be  chosen  this  way.  For  coarser 
grids  it  can  be  smaller,  which  will  yield  additional  savings,  see  [  0,  Section  4.5]. 
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Chapter  3 


Results 


3.1  Interior  problems  on  a  Cartesian  grid 

3.1.1  Schemes  of  various  accuracy 

We  first  solve  the  interior  Dirichlet  problem  for  the  constant  coefficient  homogeneous 
Helmholtz  equation: 


Ait  +  k  u  =  0  on  O, 


(3.1) 


u  |r  = 


on  the  domain  O  which  is  a  disk  of  radius  R  =  3  centered  at  the  origin.  The  auxiliary 
problem  (2.33)  is  formulated  on  a  larger  square  P  =  (—7 r,  7r)  x  (— 7r,  7r)  and  consists 
of  solving  the  inhomogeneous  Helmholtz  equation 


Av  +  k2v  =  g 
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subject  to  the  zero  Dirichlet  boundary  condition  v\a^  =  0.  To  avoid  resonances  and 
guarantee  uniqueness  of  the  solution  to  the  AP,  we  require  that  k1  2  ^  l2  +m2,  where 
l  and  m  are  any  two  integers.  The  AP  is  solved  by  a  sparse  LU  decomposition. 

We  take  the  test  solution  of  problem  (3.1)  in  the  form  of  a  plane  wave: 

u(x,  y)  =  e^kxX+kyV\  where  k2  +  k2  =  k2 ,  (3.2) 


so  that  the  boundary  data  in  (3.1)  become 


m 


_  iR(kx  cos  6 -\-ky  sin  6) 


(3.3) 


where  0  is  the  polar  angle.  The  specific  values  that  we  choose  are:  kx  =  and 
ky  =  '^k  and  k  to  be  chosen  later.  We  note  that  in  further  discussion  we  use  the 
notation  of  k  rather  than  non-dimensional  kD  where  D  is  a  radius  R  or  the  major 
elliptical  semi-axis  a  since  the  size  shape  is  fixed. 

To  discretize  the  Helmholtz  equation,  we  use  a  uniform,  in  both  directions, 
Cartesian  grid  with  size  h  on  the  square  f l  (the  domain  of  the  AP).  In  doing  so,  the 
circular  boundary  T  of  the  domain  f 2,  i.e.,  the  domain  of  the  original  problem  (3.1), 
does  not  conform  to  the  grid.  The  Helmholtz  equation  is  discretized  by  means  of 
the  following  three  schemes: 

1.  The  standard  central  difference  second  order  accurate  scheme  on  the  five- node 
stencil; 

2.  The  fourth  order  accurate  compact  scheme  of  Section  1.2.1,  which  uses  the 
nine-node  stencil  shown  in  Figure  2.2  [the  equation  in  (3.1)  is  homogeneous, 


98 


3.1.  INTERIOR  PROBLEMS  ON  A  CARTESIAN  GRID  95 

and  no  stencil  is  needed  for  the  right-hand  side].  For  the  case  of  k  =  const,  the 
scheme  simplifies  compared  to  formula  (1.8); 

3.  The  sixth  order  accurate  compact  scheme  of  [  6],  which  uses  the  same  nine- 
node  stencil. 

The  goals  of  the  computations  are  to  demonstrate  the  design  order  of  the  grid 
convergence  of  the  numerical  solution  to  the  exact  solution  for  a  non-conforming 
boundary.  We  also  determine  what  is  the  minimum  number  of  Taylor  derivatives 
needed  to  maintain  the  required  order  of  accuracy  of  the  overall  scheme;  see  discus¬ 
sion  before  formula  (2.31)  or  [  ].  An  additional  goal  is  to  show  how  the  pollution 

effect  [35,  5,  2]  manifests  itself. 

The  grid  convergence  is  studied  by  solving  on  a  sequence  of  grids  of  increasing 
dimension:  2d  x  2d ,  from  16  x  16  to  the  maximum  of  1024  x  1024.  So  for  a  given 
d  the  grid  size  is  h  =  =  ir21~d,  and  it  is  halved  every  time  the  grid  dimension  is 

increased. 

We  have  solved  problem  (3.1)  for  five  different  values  of  the  wavenumber  k:  1, 
3,  6.7,  12.8,  and  25.6.  For  the  highest  k  that  we  have  considered,  the  test  solution 
(3.2)  already  exhibits  a  fair  amount  of  oscillations  on  the  domain  fl  —  about  25  full 
wavelengths  along  the  radius,  as  shown  in  Figure  3.1.  The  results  of  computations 
for  all  k's  are  presented  in  Table  3.1  through  Table  3.5. 

Altogether,  Tables  3.1  -  3.5  show  that  for  every  scheme  we  have  tested,  the 
proposed  methodology  guarantees  the  design  rate  of  grid  convergence  for  a  non- 
conforming  boundary  and  a  Cartesian  grid. 


99 


96 


CHAPTER  3.  RESULTS 


-4  -4 

(a)  Real  Part 


-4  -4 

(b)  Imaginary  Part 


Figure  3.1:  Real  and  Imaginary  part  of  the  test  solution  (3.2)  for  k  =  25.6. 
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Grid 

Scheme 

2nd  order  ctr.  difference 

4th  order  compact 

6th  order  compact 

1 1 H  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

16  x  16 

3.474855e-2 

— 

4.590536e-4 

— 

3.626766e-6 

— 

32  x  32 

5.346252e-3 

2.7004 

5.163260e-6 

6.4742 

8.250530e-9 

8.7800 

64  x  64 

1.238241e-3 

2.1102 

1.704410e-7 

4.9209 

6.486869e-ll 

6.9908 

128  x  128 

3.001289e-4 

2.0446 

9.090205e-9 

4.2288 

1.112940e-12 

5.8651 

256  x  256 

7.389904e-5 

2.0220 

3.272063e-10 

4.7960 

2.343009e-12 

-1.0740 

512  x  512 

1.835138e-5 

2.0097 

2.457055e-ll 

3.7352 

7.287204e-12 

-1.6370 

1024  x  1024  4.571995e-6 

2.0050 

3.070920e-ll 

-0.3217 

2.411052e-ll 

-1.7262 

Table  3.1: 

Grid  convergence 

for  the  wavenumber  k 

=  1  and  the  dimension  of  the  basis 

(2.38)  M  =  17.  Note  that  the  apparent  breakdown  of  convergence  of  higher  order  schemes 
on  finer  grids  is  due  to  the  loss  of  significant  digits,  as  the  absolute  levels  of  the  error  become 
very  small  and  approach  machine  zero. 


Grid  Scheme 

2nd  order  ctr.  difference  4th  order  compact  6th  order  compact 


u  Unum||oo  rate  ||u  rtnum||oo  rate  ||u  Unum||oo  rate 


16  x  16 

2.157031 

— 

1.093211 

— 

6.252035e-2 

— 

32  x  32 

2.212195e-l 

3.2855 

1.491703e-3 

9.5174 

1.905533e-5 

11.6799 

64  x  64 

6.296501e-2 

1.8129 

4.695925e-5 

4.9894 

2.743013e-7 

6.1183 

128  x  128 

1.621645e-2 

1.9571 

2.736886e-6 

4.1008 

3.956555e-9 

6.1154 

256  x  256 

4.049416e-3 

2.0017 

1.612331e-7 

4.0853 

5.830238e-ll 

6.0845 

512  x  512 

1.008930e-3 

2.0049 

9.823236e-9 

4.0368 

1.288003e-12 

5.5003 

1024  x  1024  2.515190e-4 

2.0041 

6.235303e-10 

3.9777 

7.870095e-12 

-2.6112 

Table  3.2: 

Grid  convergence 

for  the  wavenumber  k 

=  3  and  the  dimension  of  the  basis 

(2.38)  M  =  28. 


The  dimension  M  of  the  basis  (2.38)  is  chosen  by  Fourier  transforming  the 
boundary  data  (3.3)  of  problem  (3.1)  and  truncating  the  series  at  the  machine 
precision  level  (double  precision).  The  resulting  values  of  M  for  every  k  are 
provided  in  the  captions  to  Tables  3.1  -  3.5.  We  see  that  M  increases  as  k  increases. 
This  is  not  surprising,  as  the  solution  becomes  more  oscillatory;1  for  example,  our 
highest  k  =  25.6  corresponds  to  over  75  full  wavelengths  along  the  circumference 

1  Convergence  of  the  Fourier  series  remains  exponential  due  to  the  smoothness,  but  the  constants 
become  larger. 
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Grid  Scheme 

2nd  order  ctr.  difference  4th  order  compact  6th  order  compact 


u  tinum||oo  rate  ll'U  Unum||oo  rate  ||w  Unum||oo  rate 


16  x  16 

8.228113 

— 

7.976387 

— 

1.459757e+l 

— 

32  x  32 

4.299483 

0.9364 

1.001469e-l 

6.3155 

1.071898e-2 

10.4113 

64  x  64 

1.933134 

1.1532 

4.111705e-3 

4.6062 

1.345917e-4 

6.3154 

128  x  128 

3.065574e-2 

2.6567 

2.420283e-4 

4.0865 

1.845635e-6 

6.1883 

256  x  256 

7.028536e-3 

2.1249 

1.488596e-5 

4.0232 

2.757929e-8 

6.0644 

512  x  512 

1.861682e-3 

1.9166 

9.101549e-7 

4.0317 

4.192718e-10 

6.0396 

1024  x  1024  4.726287e-4 

1.9778 

5.640010e-8 

4.0123 

1.170244e-ll 

5.1630 

Table  3.3: 

Grid  convergence 

for  the  wavenumber  k 

=  6.7  and  the  dimension  of  the  basis 

(2.38)  M  =  43. 


Grid  Scheme 


2nd  order  ctr.  difference  4th  order  compact  6th  order  compact 


1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

16  x  16 

3.135403e+l 

— 

7.284196e+l 

— 

9.174488e+l 

— 

32  x  32 

2.693366e+l 

0.2192 

4.960223 

3.8763 

1.344958e+l 

2.7701 

64  x  64 

8.177246 

0.2192 

1.233802 

2.0073 

8.032610e-2 

7.3875 

128  x  128 

1.095035e+l 

-0.4213 

3.200884e-2 

5.2685 

1.039313e-3 

6.2722 

256  x  256 

2.603452 

2.0725 

2.048553e-3 

3.9658 

1.395774e-5 

6.2184 

512  x  512 

6.781712e-l 

1.9407 

1.277844e-4 

4.0028 

2.125559e-7 

6.0371 

1024  x  1024 

1.448771e-l 

2.2268 

7.718401e-6 

4.0493 

3.172309e-9 

6.0662 

Table  3.4:  Grid 

convergence 

for  the  wavenumber  k  = 

12.8  and  the  dimension  of  the  basis 

(2.38)  M  =  66. 


R  =  3.  On  the  other  hand,  we  also  see  in  Tables  3.1  through  3.5  that  the  accuracy 
actually  achieved  on  the  grid  is  often  orders  of  magnitude  less  than  the  machine 
precision.  This  indicates  that  the  chosen  M  may  be  superfluous,  and  the  same 
accuracy  of  the  solution  can  be  obtained  using  a  smaller  basis  (2.38)  at  a  lower 
computational  cost.  For  example,  the  fourth  and  sixth  order  computations  presented 
in  Table  3.5  (k  =  25.6)  were  repeated  for  various  values  of  M  with  similar  results. 
From  Tables  3.6  and  3.7  one  learns  that  for  the  finest  grids  the  higher  accuracy 
reached  with  the  large  value  of  M  degrades  when  the  value  of  M  decreased  more 
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then  20  percent.  This  result  is  not  surprising  since  Fourier  series  were  cut  at  10~2 
and  10-3  for  M  =  90  and  M  =  85  respectively. 


Grid 

Scheme 

2nd  order  ctr. 

difference 

4th  order  compact 

6th  order  compact 

||^  Unum||oo 

rate 

1 1  ^  ^num  1 1  oo 

rate 

1 1  ^  Unum  1 1  oo 

rate 

16  x  16 

1.144291e+2 

— 

2.065713e+3 

— 

1.175546e+4 

— 

32  x  32 

4.851901e+l 

1.2378 

8.885777e+l 

4.5390 

4.453015e+l 

8.0443 

64  x  64 

1.280298e+l 

1.9221 

2.818431 

4.9785 

1.219721e+l 

1.8682 

128  x  128 

1.901798e+l 

-0.5709 

4.128656e-l 

2.7711 

1.039313e-3 

5.5781 

256  x  256 

1.448009e+l 

0.3933 

1.737760e-l 

1.2484 

1.973019e-3 

7.0158 

512  x  512 

4.563927 

1.6657 

4.317500e-3 

5.3309 

2.883989e-5 

6.0962 

1024  x  1024 

3.892365 

0.2296 

2.603055e-4 

4.0519 

4.398634e-7 

6.0349 

Table  3.5:  Grid  convergence  for  the  wavenumber  A;  =  25.6  and  the  dimension  of  the  basis 
(2.38)  M  =  111. 


Grid 

M 

111 

100 

90 

85 

11^  Unum  1 1  oo 

1 1  ^  Unum  1 1  oo 

1 1  ^  Unum  1 1  oo 

1 1  ^  Unum  1 1  oo 

16  x  16 

2.065713e+3 

2.253702e+3 

2.222923e+3 

2.222930e+3 

32  x  32 

8.885777e+l 

6.484520e+l 

6.827580e+l 

5.631450e+l 

64  x  64 

2.818431 

2.813891 

2.503964 

2.192968e+0 

128  x  128 

4.128656e-l 

4.070042e-l 

4.048406e-l 

4.057072e-l 

256  x  256 

1.737760e-l 

1.735016e-l 

1.733309e-l 

1.734153e-l 

512  x  512 

4.317500e-3 

4.317889e-3 

4.317991e-3 

1.220778e-2 

1024  x  1024 

2.603055e-4 

2.603074e-4 

5.893227e-4 

1.066776e-2 

Table  3.6:  Behavior  of  the  schemes  for  various  dimensions  of  the  basis  (2.38)  M  -  4th 
order  compact  scheme.  The  wavenumber  is  A  =  25.6. 


Grid  M 

111  100  90  85 


Unum  1 1  oo  1 1  ^  Unum  1 1  oo  1 1 U  Unum  1 1  oo  1 1 U  Unum  1 1  oo 


16  x  16 

1.175546e+4 

1.352850e+4 

1.276175e+4 

1.237621e+4 

32  x  32 

4.453015e+l 

4.953919e+l 

4.407588e+l 

4.927018e+l 

64  x  64 

1.219721e+l 

7.862976 

7.118779 

6.783102 

128  x  128 

1.039313e-3 

2.549135e-l 

2.530015e-l 

2.529428e-l 

256  x  256 

1.973019e-3 

1.973351e-3 

1.971636e-3 

1.491697e-2 

512  x  512 

2.883989e-5 

2.884148e-5 

6.756577e-4 

1.194734e-2 

1024  x  1024 

4.398634e-7 

4.882227e-7 

5.914269e-4 

1.065307e-2 

Table  3.7:  Behavior  of  the  schemes  for  various  dimensions  of  the  basis  (2.38)  M  -  6th 
order  compact  scheme.  The  wavenumber  is  A  =  25.6. 
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Grid 

Scheme 

2nd  order  ctr. 

difference 

4th  order  compact 

6th  order  compact 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

16  x  16 

3.752407 

— 

2.907471e-l 

— 

7.342967e-2 

— 

32  x  32 

2.290364e-l 

4.0342 

1.010567e-2 

4.8465 

1.432143e-4 

9.0020 

64  x  64 

6.503502e-2 

1.8163 

8.516129e-4 

3.5688 

3.976092e-6 

5.1707 

128  x  128 

1.990979e-2 

1.7077 

7.105160e-5 

3.5833 

9.373316e-8 

5.4066 

256  x  256 

5.743664e-3 

1.7934 

4.898880e-6 

3.8583 

1.693751e-9 

5.7903 

512  x  512 

2.130201e-3 

1.4310 

3.736573e-7 

3.7127 

3.099734e-ll 

5.7719 

1024  x  1024 

4.835611e-4 

2.1392 

2.310974e-8 

4.0151 

7.720891e-12 

2.0053 

Table  3.8:  The  deterioration  in  grid  convergence  for  the  wavenumber  k  =  3  and  the  dimen¬ 
sion  of  the  basis  (2.38)  M  =  28  when  n  —  1  degree  Taylor  extension  used. 


Previously,  in  Section  2. 3. 1.2  we  claimed  that  in  order  to  obtain  nth  order  ap¬ 
proximation  to  the  original  problem  one  needs  Equation  Based  Taylor  extension  of 
nth  degree.  Thus,  we  note  that  all  the  computations  presented  in  Tables  3.1  -  3.5 
and  in  Tables  3.6,  3.7  were  conducted  using  a  nth  degree  Equation  Based  Taylor 
extension.  Theoretically,  the  degree  of  Taylor  extension  were  predicted  to  be  n  +  2, 
see  [  ].  To  see  whether  or  not  our  current  (lower)  choice  of  Taylor’s  order  can  be 

improved  further,  we  have  conducted  similar  computations,  but  for  an  even  smaller, 
n  —  1,  degree  of  Taylor  extension.  In  Table  3.8,  we  present  the  results  for  k=  3. 
The  data  show  a  certain  deterioration  of  the  convergence  rate  (cf.  Table  3.2),  which 
indicates  that  the  number  of  terms  in  the  equation-based  extension  formulae  should 
not  be  taken  any  lower  than  the  order  of  accuracy  n  of  the  scheme  (i.e.  nth  degree 
Taylor  extension). 
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3.1.2  Variable  wavenumber  Helmholtz  equation  with  fourth  order 
accuracy 

We  now  use  the  fourth  order  accurate  compact  scheme  (1.8)  to  solve  the  inhomoge¬ 
neous  Helmholtz  equation  (1.1)  with  a  variable  wavenumber  inside  the  circles  and 
ellipses,  subject  to  Dirichlet  or  Neumann  boundary  conditions.  The  goal  of  the 
computations  is  to  demonstrate  the  capability  of  the  proposed  method  to  address 
variable  coefficients  and  various  types  of  the  boundary  conditions,  and  again,  to 
show  the  design  order  of  grid  convergence  for  non-conforming  boundaries  (the  dis¬ 
cretization  grid  is  always  Cartesian).  The  domain  fh  is  either  a  disk  of  radius  R  =  1 
centered  at  the  origin,  or  the  interior  of  the  ellipse  with  the  major  semi-axis  a  =  1 
and  minor  semi-axis  b  =  1/2,  see  formula  (1.5). 

The  Helmholtz  equation  (1.1)  that  we  solve  on  the  domain  fii  has  a  variable 
wavenumber  k.  For  the  case  of  the  disk  we  choose 

k  =  kQe-1^r~r  °)6r6c08*,  (3.4) 

and  for  the  case  of  the  ellipse  we  take 

k  =  koe-10^-^6 ,  (3.5) 

where  r  is  the  polar  radius  and  9  is  the  polar  angle  and  the  parameter  ?’o  =  1.6.  The 
profiles  of  k  are  schematically  shown  in  Figure  3.2. 
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(b)  fli  is  an  ellipse,  formula  (  3.5) 


Figure  3.2:  Profiles  of  the  variable  wavenumber  k  on  for  =  25;  the  part  inside  fli  is 
emphasized. 


In  either  case,  circle  or  ellipse,  the  exact  solution  is  chosen  in  the  form: 


u  = 


^ikx 


(3.6) 
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Since  k  is  variable,  see  formulae  (3.4)  and  (3.5),  this  solution  is  not  a  plane  wave,  as 
shown  in  Figure  3.3  and  3.4.  The  corresponding  right-hand  side  f(x,y )  in  formula 
(1.1)  is  obtained  by  backward  engineering,  i.e.,  by  substituting  u  given  by  (3.6)  into 
the  left-hand  side  of  the  Helmholtz  equation. 

The  boundary  condition  at  T  =  for  the  Helmholtz  equation  (1.1)  can  be  of 
either  Dirichlet  or  Neumann  type.  The  required  boundary  data  are  also  obtained 
by  backward  engineering,  i.e.,  by  taking  the  trace  of  either  the  solution  u  itself  or 
its  normal  derivative  at  the  boundary  T. 

When  Hi  is  a  disk  of  radius  R  =  1,  the  AP  (see  2.33)  is  formulated  on  the  square 

=  {(-T  y)  I  -  1.2  <  X,  y  ^  1.2} 

with  the  following  boundary  conditions:  v  =  0  at  y  =  ±1.2, and 

dv  ,  dv  , 

— — b  iv  =  0  at  x  =  1.2  and  — - iv  =  0  at  x  =  —1.2.  (3-7) 

dx  dx 

The  pair  of  complex  boundary  conditions  (3.7)  guarantees  that  regardless  of  k  there 
will  be  no  resonances  in  the  solution  of  the  AP  on  the  square  Ho- 

In  the  case  of  the  ellipse,  the  AP  is  formulated  on  the  rectangle 
H0  =  {{x,y)  |  -  1.2  <  s  <  1.2,  -0.7  <  y  ^  0.7} 
with  the  boundary  conditions: 


v  =  0  at  y  =  ±0.7, 
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-1.5  -1.5 

(a)  Real  part 


-1.5  -1.5 

(b)  Imaginary  part 


Figure  3.3:  Real  and  Imaginary  part  of  the  test  solution  (3.6)  in  circle  for  fco  =  25. 


and  the  same  complex  boundary  conditions  (3.7)  at  x  =  ±1.2. 


Similar  to  Section  3.1.1,  the  AP  is  discretized  on  a  sequence  of  uniform,  in  each 
direction,  Cartesian  grids  of  dimension  2d  x  2rf,  with  a  maximum  of  2048  x  2048. 
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(a)  Real  Part 


-1 

■1  -1.5 

(b)  Imaginary  Part 


Figure  3.4:  Real  and  Imaginary  part  of  the  test  solution  (3.6)  in  ellipse  for  k0  =  25. 


For  a  given  d,  the  grid  size  in  the  case  of  a  square  is  h  =  ,  and  the  grid  sizes  in 

the  case  of  a  rectangle  are  hx  =  and  hy  =  \^.  The  grid  sizes  are  halved  every 
time  d  is  incremented  by  1,  which  is  convenient  for  studying  the  convergence.  As  in 


109 


106 


CHAPTER  3.  RESULTS 


Section  3.1.1,  the  AP  is  also  solved  by  the  sparse  LU  decomposition. 

Numerical  results  for  solving  the  Dirichlet  problem  for  the  variable  coefficient 
Helmholtz  equation  (1.1)  are  presented  in  Table  3.9  in  the  case  of  a  circle  and  in 
Table  3.10  in  the  case  of  an  ellipse  of  an  aspect  ratio  AR  =  2  for  the  range  of  fc’s 
that  we  have  investigated,  ko=  5,  15,  and  25.  In  Table  3.11  we  fix  the  wavenumber 
to  be  ko=  10  and  vary  the  aspect  ratio  of  an  ellipse  for  the  values  AR  =  2,  4  and 
8.  Ellipses  with  even  higher  AR  are  presented  in  later  sections.  We  stress  that  the 
data  in  the  tables  fully  corroborate  the  design  fourth  order  rate  of  grid  convergence 
for  the  compact  Cartesian  scheme  (1.8)  when  the  non-conforming  boundaries  are 
handled  by  the  method  of  difference  potentials. 


Similar  numerical  results  for  the  Neumann  problem  are  presented  in  Table  3.12 
for  the  circle,  in  Table  3.13  for  the  ellipse,  and  in  Table  3.14  for  ellipses  with  different 
aspect  ratios.  As  with  the  Dirichlet  problem,  the  data  in  the  tables  fully  corroborate 
the  design  fourth  order  rate  of  grid  convergence  of  the  proposed  methodology. 


Grid 

Circle 

k0  =  5,  M 

=  42 

k0  =  15,  M 

=  54 

k0  =  25,  M 

=  67 

1 1  ^  W-num  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

16  x  16 

3.038830 

— 

1.517016e+l 

— 

1.571341e+2 

— 

32  x  32 

2.231693e-2 

7.0892 

3.324707 

5.5119 

1.304138e+l 

3.5908 

64  x  64 

1.405180e-3 

3.9893 

1.032080e-l 

5.0096 

2.744285 

2.2486 

128  x  128 

7.302520e-5 

4.2662 

5.746378e-3 

4.1668 

5.751020e-2 

5.5765 

256  x  256 

4.465171e-6 

4.0316 

3.454849e-4 

4.0560 

3.678247e-3 

3.9667 

512  x  512 

2.701632e-7 

4.0468 

2.125100e-5 

4.0230 

2.265488e-4 

4.0211 

1024  x  1024 

1.680068e-8 

4.0072 

1.321587e-6 

4.0072 

1.405292e-5 

4.0109 

2048  x  2048 

1.040726e-9 

4.0129 

8.239623e-8 

4.0035 

8.745530e-7 

4.0062 

Table  3.9:  Grid  convergence  of  the  solution  to  the  Dirichlet  problem  for  the  circle  R=  1. 
Variable  coefficient  Helmholtz  equation  (1.1)  and  a  fourth  order  compact  scheme  (1.8). 
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Grid 

Ellipse 

k0  =  5,  M 

=  42 

k0  =  15,  M 

=  54 

k0  =  25,  M 

=  67 

1 1  ^  M-num  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  l^num  1 1  oo 

rate 

16  x  16 

8.283939 

— 

3.950658e+l 

— 

5.206655e+l 

— 

32  x  32 

1.978583e-2 

8.7097 

6.015599e-l 

6.0372 

3.388063e+l 

0.6199 

64  x  64 

3.104902e-4 

5.9938 

7.001777e-3 

6.4248 

8.662025e-2 

8.6115 

128  x  128 

1.659692e-5 

4.2256 

7.492233e-4 

3.2243 

5.811711e-3 

3.8977 

256  x  256 

5.597237e-7 

4.8901 

2.551093e-5 

4.8762 

3.104959e-4 

4.2263 

512  x  512 

2.094249e-8 

4.7402 

1.551669e-6 

4.0392 

1.881038e-5 

4.0450 

1024  x  1024 

6.565249e-10 

4.9954 

9.538440e-8 

4.0239 

1.160326e-6 

4.0189 

2048  x  2048 

2.761463e-ll 

4.5713 

5.897927e-9 

4.0155 

7.200706e-8 

4.0102 

Table  3.10:  Grid  convergence  of  the  solution  to  the  Dirichlet  problem  for  the  ellipse  o  =  l, 
b=  I.  Variable  coefficient  Helmholtz  equation  (1.1)  and  a  fourth  order  compact  scheme 
(1.8). 


Grid 

M  =  54,  AR 

=  2 

M  =  72,  AR  = 

4 

M  =  98,  AR  = 

8 

||uh-u2,l||00 

rate 

||uh-u2h||00 

rate 

HtV-u^lU 

rate 

16  x  16 

3.707995e+l 

— 

6.563552  e+2 

— 

3.514365e+2 

— 

32  x  32 

1.678518e-l 

7.7873 

3.475799  e+2 

0.9171 

4.113777e+4 

-6.8711 

64  x  64 

1.771571e-3 

6.5660 

2.900153  e+3 

-3.0607 

1.338930e+l 

1.5852 

128  x  128 

1.915687e-4 

3.2091 

2.715441  e-4 

23.3484 

1.148908e+2 

-3.1011 

256  x  256 

5.270431e-6 

5.1838 

1.038533  e-5 

4.7086 

2.092065e-3 

15.7450 

512  x  512 

2.157025e-7 

4.6108 

8.242304  e-7 

3.6554 

1.650560e-6 

10.3078 

1024  x  1024 

1.315849e-8 

4.0350 

2.520539  e-8 

5.0312 

7.926636e-8 

4.3801 

2048  x  2048 

7.490176e-10 

4.1349 

6.660874  e-10 

5.2419 

6.770679e-9 

3.5493 

Table  3.11:  Grid  convergence  of  the  solution  to  the  Dirichlet  problem  for  the  wavenumber 
fc  =  10  and  the  ellipses  a=l,  b£  Variable  coefficient  Helmholtz  equation  (1.1)  and 

a  fourth  order  compact  scheme  (1.8). 


Grid 

Ellipse 

k0  =  5,  M 

=  42 

k0  =  15,  M 

=  54 

k0  =  25,  M 

=  67 

1 1  ^  Unum  1 1  oo 

rate 

1 1 U  Unum  1 1  oo 

rate 

1 1 V*  l^num  1 1  oo 

rate 

16  x  16 

4.267776e+l 

— 

1.772119e+2 

— 

3.649132e+2 

— 

32  x  32 

9.712126e-2 

8.7795 

2.316084 

6.2576 

3.309569e+l 

3.4628 

64  x  64 

7.548180e-3 

3.6856 

6.203786e-2 

5.2224 

1.573478e-2 

7.7165 

128  x  128 

4.486249e-4 

4.0725 

4.713176e-3 

3.7184 

1.589413e-3 

3.3074 

256  x  256 

2.486193e-5 

4.1735 

2.419222e-4 

4.2841 

6.383346e-4 

4.6380 

512  x  512 

1.372890e-6 

4.1787 

1.635393e-5 

3.8868 

4.329843e-5 

3.8819 

1024  x  1024 

9.028545e-8 

3.9266 

9.750050e-7 

4.0681 

2.425497e-6 

4.1580 

2048  x  2048 

5.198146e-9 

4.1184 

6.308512e-8 

3.9500 

1.584606e-8 

3.9361 

Table  3.13:  Grid  convergence  of  the  solution  to  the  Neumann  problem  for  the  ellipse  a=  1, 
b=  I.  Variable  coefficient  Helmholtz  equation  (1.1)  and  a  fourth  order  compact  scheme 
(1.8). 
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Grid 

Circle 

k0  =  5,  M 

=  42 

k0  =  15,  M 

=  54 

k0  =  25,  M 

=  67 

1 1 U  Unum  1 1  oo 

rate 

1 1 U  Itnum  1 1  oo 

rate 

1 1  Unum  1 1  oo 

rate 

16  x  16 

1.188942 

— 

8.839874e+l 

— 

1.713283e+2 

— 

32  x  32 

1.801846e-2 

6.0441 

6.391570e+l 

5.5119 

7.346974e+l 

1.2215 

64  x  64 

1.245872e-3 

3.8542 

3.179250e-l 

5.0096 

2.938586 

4.6440 

128  x  128 

5.731111e-5 

4.4422 

3.336937e-2 

4.1668 

4.978325e-2 

5.8833 

256  x  256 

4.343596e-6 

3.7219 

2.037841e-3 

4.0560 

2.970674e-3 

4.0668 

512  x  512 

2.118921e-7 

4.3575 

1.218354e-4 

4.0230 

1.810039e-4 

4.0367 

1024  x  1024 

1.467516e-8 

3.8519 

7.567901e-6 

4.0072 

1.144116e-5 

3.9837 

2048  x  2048 

8.996365e-10 

4.0279 

4.647261e-7 

4.0035 

7.144455e-7 

4.0013 

Table  3.12:  Grid  convergence  of  the  solution  to  the  Neumann  problem  for  the  circle  R=  1. 
Variable  coefficient  Helmholtz  equation  (1.1)  and  a  fourth  order  compact  scheme  (1.8). 


Grid 

M  =  54,  AR  = 

2 

M  =  72,  AR  = 

4 

M  =  98,  AR  = 

8 

||u^  -  u2h||00 

rate 

llu^-U^Hoc 

rate 

!|uh'-u2'l||00 

rate 

16  x  16 

1.118772e+2 

— 

5.034143e+3 

— 

2.361988e+4 

— 

32  x  32 

4.950664e-l 

7.8201 

1.350765e+4 

-1.4240 

2.316460e+5 

-3.2938 

64  x  64 

9.400937e-3 

5.7187 

6.645722e+l 

7.6671 

1.798778e+l 

13.6526 

128  x  128 

7.042638e-4 

3.7386 

2.195422e-3 

14.8856 

5.980095e+l 

-1.7332 

256  x  256 

3.226175e-5 

4.4482 

1.562795e-4 

3.8123 

3.296630e-2 

10.8250 

512  x  512 

2.080956e-6 

3.9545 

1.028099e-5 

3.9261 

4.720918e-5 

9.4477 

1024  x  1024 

1.202756e-7 

4.1128 

5.577696e-7 

4.2042 

3.309864e-6 

3.8342 

2048  x  2048 

7.788341e-9 

3.9489 

3.868259e-8 

3.8499 

2.209241e-7 

3.9051 

Table  3.14:  Grid  convergence  of  the  solution  to  the  Neumann  problem  for  the  wavenumber 
/c  =  10  and  the  ellipses  o=l,  &€  Variable  coefficient  Helmholtz  equation  (1.1)  and 

a  fourth  order  compact  scheme  (1.8). 


3.2  Exterior  scattering  problems 

We  next  consider  the  scattering  of  an  incoming  plane  wave  with  a  given  frequency 
(wavelength)  and  given  angle  of  incidence  off  an  elliptical  body  with  a  given  aspect 
ratio.  In  our  simulations,  we  take  the  major  semi-axis  of  the  ellipse  to  be  a  =  1.8, 
while  its  minor  semi-axis  varies  between  b  =  0.9  and  b  =  0.18,  which  yields  aspect 
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ratios  between  2  and  10.  The  wavenumber  in  the  Helmholtz  equation  (  1.2a)  varies 
between  =  1  and  ko  =  25,  which  yields  the  variation  of  the  wavelength  between 
Ao  =  2-7T  and  Ao  =  27r/25,  i.e.,  between  roughly  twice  the  size  2a  of  the  ellipse  and 
about  8%  of  this  size.  We  consider  several  values  of  the  angle  of  incidence  between 
0°  and  50°  with  respect  to  the  direction  of  the  major  axis.  We  also  consider  both 
Dirichlet  and  Neumann  boundary  conditions  at  the  contour  T,  i.e.,  at  the  perimeter 
of  the  ellipse.  In  the  context  of  acoustics,  the  former  corresponds  to  sound-soft 
scattering,  whereas  the  latter  corresponds  to  sound-hard  scattering.  The  exterior 
AP  is  solved  on  the  domain  Hi  shaped  as  an  annulus,  Hi  =  {Rq  ^  r  ^  i?i},  see 
Figure  2.4,  with  Rq  that  may  vary  between  0.1  (for  larger  aspect  ratios)  and  0.3 
(for  smaller  aspect  ratios),  and  R±  =  2.  This  AP  is  discretized  on  a  uniform,  in 
each  direction,  polar  grid  that  may  have  between  64  x  64  and  4096  x  4096  cells. 
The  quantity  M  that  represents  the  dimension  of  the  basis  on  T,  see  formula  (2.38), 
is  grid-independent  and  chosen  so  as  to  guarantee  that  the  given  boundary  data 
(Dirichlet  or  Neumann)  are  approximated  by  the  corresponding  finite  Fourier  series 
up  to  the  machine  precision.  The  problem  is  solved  using  the  simplified  methodology 
of  Section  2.3.6.  In  doing  so,  the  discrete  exterior  AP  is  integrated  by  means  of  the 
separation  of  variables  combined  with  a  FFT.  The  exact  nonlocal  ABC  at  the  outer 
circle  r  =  R±  is  conveniently  set  in  the  Fourier  space,  see  [  ]. 

As  the  overall  set  of  results  for  all  wavenumbers,  incidence  angles,  aspect  ratios, 
etc.,  is  rather  large,  we  have  chosen  to  show  only  a  representative  sample.  In  Fig¬ 
ures  3.5  and  3.6,  we  show  the  schematic  geometry  for  two  ellipses  —  of  aspect  ratio 
2  and  of  aspect  ratio  10  (cf.  Figure  2.4). 


To  assess  the  grid  convergence,  we  do  not  evaluate  the  exact  solution  using  its 
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Figure  3.5:  Schematic  of  the  polar  grid  for  the  exterior  AP,  the  elliptic  scatterer  of  aspect 
ratio  2,  and  the  grid  boundary  7. 

expansion  with  respect  to  Mathieu  functions  [6],  because  this  may  entail  numerical 
difficulties  of  its  own.  We  rather  evaluate  the  infinity  norm  of  the  difference  between 
the  numerical  solutions  obtained  on  two  consecutive  grids,  uh  and  u2h . 

Tables  3.15  through  3.22  demonstrate  the  design  fourth  order  rate  of  the  grid 
convergence  for  the  case  of  a  Dirichlet  boundary  condition  on  T.  We  note  that  the 
convergence  on  coarser  grids  looks  somewhat  more  “erratic”  for  slenderer  ellipses. 
This  is  likely  accounted  for  by  insufficient  grid  resolution  in  the  areas  of  high  curva¬ 
ture,  i.e.,  near  the  tips  of  the  major  axis.  Nonetheless,  on  finer  grids  the  convergence 
rate  approaches  its  correct  asymptotic  value  of  4.  Similar  results  are  obtained  for 
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Figure  3.6:  Schematic  of  the  polar  grid  for  the  exterior  AP,  the  elliptic  scatterer  of  aspect 
ratio  10,  and  the  grid  boundary  7. 

the  Neumann  boundary  condition  on  T,  see  Tables  3.24  and  3.25,  as  well  as  for  a 
variety  of  other  Dirichlet  and  Neumann  cases  that  are  not  presented  in  these  tables. 


Grid 

CM 

II 

1—1 

II 

O 

-se 

k0  =  10,  M  = 

37 

k0  =  25,  M  = 

69 

|uh  -  U^lloo 

rate 

Hu'1  -u2h||00 

rate 

|uh  —  u2,l||00 

rate 

64  x  64 

2.737716ee-3 

— 

6.104800 

— 

3.945149e+l 

— 

128  x  128 

3.658040ee-4 

2.9038 

6.761083e-2 

6.4965 

1.020563e+2 

-1.3712 

256  x  256 

2.075372ee-5 

4.1396 

3.888347e-3 

4.1200 

3.843873e-l 

8.0526 

512  x  512 

1.134933ee-6 

4.1927 

2.382304e-4 

4.0287 

2.189631e-2 

4.1338 

1024  x  1024 

6.890566ee-8 

4.0418 

1.481515e-5 

4.0072 

1.329834e-3 

4.0414 

2048  x  2048 

4.235718ee-9 

4.0239 

9.247807e-7 

4.0018 

8.251 120e-5 

4.0105 

4096  x  4096 

2.637849ee-10 

4.0052 

5.778413e-8 

4.0004 

5.147419e-6 

4.0027 

Table  3.15:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  0°  about  an  ellipse 
with  aspect  ratio  2. 
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Grid 

k0  =  1,  M  = 

12 

k0  =  10,  M  = 

37 

k0  =  25,  M  = 

69 

||u',-U2h||0o 

rate 

||uh-u2,l||00 

rate 

||uh-U2h||oo 

rate 

64  x  64 

2.924441e-3 

— 

1.399531 

— 

1.588919e+l 

— 

128  x  128 

3.686772e-4 

2.9877 

9.024242e-2 

3.9550 

2.571935e+l 

-0.6948 

256  x  256 

2.092322e-5 

4.1392 

5.042251e-3 

4.1617 

4.528627e-l 

5.8276 

512  x  512 

1.140182e-6 

4.1978 

3.077507e-4 

4.0342 

2.614107e-2 

4.1147 

1024  x  1024 

7.045679e-8 

4.0164 

1.912254e-5 

4.0084 

1.576210e-3 

4.0518 

2048  x  2048 

4.352290e-9 

4.0169 

1.193405e-6 

4.0021 

9.766250e-5 

4.0125 

4096  x  4096 

2.708936e-10 

4.0060 

7.456103e-8 

4.0005 

6.090803e-6 

4.0031 

Table  3.16:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  35°  about  an  ellipse 
with  aspect  ratio  2. 


Grid 

ko  =  1,  M  = 

12 

k0  =  10,  M  = 

37 

k0  =  25,  M  = 

69 

||uh-U2h||oo 

rate 

Uh  —  U2'*  oo 

rate 

||uh-U2h||oo 

rate 

64  x  64 

3.647013e-2 

— 

6.942957 

— 

6.094755e+5 

— 

128  x  128 

3.229204e-3 

3.4975 

1.556797e-l 

5.4789 

8.271102 

16.1691 

256  x  256 

5.948035e-4 

2.4407 

5.275428e-3 

4.8831 

4.000126e-l 

4.3700 

512  x  512 

2.733103e-5 

4.4438 

2.869505e-4 

4.2004 

2.340879e-2 

4.0949 

1024  x  1024 

1.429272e-6 

4.2572 

1.782570e-5 

4.0088 

1.410523e-3 

4.0527 

2048  x  2048 

8.252679e-8 

4.1143 

1.112633e-6 

4.0019 

8.738691e-5 

4.0127 

4096  x  4096 

5.198425e-9 

3.9887 

6.951841e-8 

4.0004 

5.449762e-6 

4.0032 

Table  3.17:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  35°  about  an  ellipse 
with  aspect  ratio  3. 


Grid 

II 

II 

O 

12 

k0  =  10,  M  = 

37 

k0  =  25,  M  = 

69 

|llh  —  U^llco 

rate 

||u'1  -  u2'l||00 

rate 

|  -  U2^1 1 1  oo 

rate 

64  x  64 

2.261699e-2 

— 

5.730228 

— 

1.742757e+3 

— 

128  x  128 

3.149537e-3 

2.8442 

9.782407e-2 

5.8723 

9.548658 

7.5119 

256  x  256 

6.271734e-4 

2.3282 

4.605138e-3 

4.4089 

2.276872 

2.0682 

512  x  512 

2.917059e-5 

4.4263 

2.738948e-4 

4.0716 

2.445093e-2 

6.5410 

1024  x  1024 

1.537777e-6 

4.2456 

1.701433e-5 

4.0088 

1.483170e-3 

4.0431 

2048  x  2048 

8.870759e-8 

4.1156 

1.061794e-6 

4.0022 

9.201076e-5 

4.0107 

4096  x  4096 

5.523241e-9 

4.0055 

6.633567e-8 

4.0006 

5.739903e-6 

4.0027 

Table  3.18:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  15°  about  an  ellipse 
with  aspect  ratio  3. 
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Grid 

k0  =  1,  M  = 

13 

k0  =  10,  M  = 

39 

k0  =  25,  M  = 

73 

IK  -  u2'l||00 

rate 

IK  -  U^Hoo 

rate 

I  K  -  U2h||oo 

rate 

64  x  64 

7.024091e+l 

— 

1.800415e+3 

— 

1.189542e+7 

— 

128  x  128 

1.117795e+l 

2.6517 

1.080259e+l 

7.3808 

2.905446e+3 

11.9994 

256  x  256 

8.069814e-3 

10.4358 

5.008329e-2 

7.7528 

6.346536e-l 

12.1605 

512  x  512 

1.523137e-3 

2.4055 

5.161442e-3 

3.2785 

2.456628e-2 

4.6912 

1024  x  1024 

7.604331e-5 

4.3241 

4.003719e-4 

3.6884 

1.471074e-3 

4.0617 

2048  x  2048 

3.763327e-6 

4.3367 

1.942630e-5 

4.3653 

9.124847e-5 

4.0109 

4096  x  4096 

2.072289e-7 

4.1827 

1.066348e-6 

4.1873 

5.691707e-6 

4.0029 

Table  3.19:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  15°  about  an  ellipse 
with  aspect  ratio  5. 


Grid 

ko  =  1,  M  = 

13 

k0  =  10,  M  = 

39 

k0  =  25,  M  = 

73 

K  -  co 

rate 

IK-U^llco 

rate 

IK-U^Hoo 

rate 

64  x  64 

7.754767e+l 

— 

1.618151e+3 

— 

2.224965e+4 

— 

128  x  128 

1.489596e+l 

2.3802 

4.674100 

8.4354 

2.185893e+l 

9.9913 

256  x  256 

1.111821e-2 

10.3878 

1.362536e-l 

5.1003 

7.391991e-l 

4.8861 

512  x  512 

1.435914e-3 

2.9529 

5.846135e-3 

4.5427 

5.567719e-2 

3.7308 

1024  x  1024 

6.645732e-5 

4.4334 

4.166078e-4 

3.8107 

1.532545e-3 

5.1831 

2048  x  2048 

3.288638e-6 

4.3369 

2.038867e-5 

4.3529 

9.506068e-5 

4.0109 

4096  x  4096 

1.804245e-7 

4.1880 

1.145084e-6 

4.1542 

5.933124e-6 

4.0020 

Table  3.20:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  50°  about  an  ellipse 
with  aspect  ratio  5. 


Grid 

k0  =  1,  M  = 

11 

k0  =  10,  M  = 

32 

k0  =  25,  M  = 

56 

lK-U2,l||co 

rate 

iK-U^Hoo 

rate 

lK-U2h||oo 

rate 

64  x  64 

2.247929e+2 

— 

4.633311e+3 

— 

1.553633e+5 

— 

128  x  128 

5.212738e+2 

-1.2134 

4.70154e+2 

3.3008 

1.547453e+5 

0.0057 

256  x  256 

8.031341e+2 

-0.6236 

4.419326e+2 

0.0893 

8.284062e+3 

7.5453 

512  x  512 

1.195681e-2 

16.0355 

4.052018e-l 

10.0910 

9.340049e-l 

9.7927 

1024  x  1024 

4.655482e-3 

1.3608 

2.785232e-2 

3.8628 

8.518332e-2 

3.4548 

2048  x  2048 

5.918121e-4 

2.9757 

1.895585e-3 

3.8771 

2.569198e-3 

5.0512 

4096  x  4096 

2.142775e-5 

4.7876 

8.621134e-5 

4.4586 

1.937799e-4 

3.7288 

Table  3.21:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  50°  about  an  ellipse 
with  aspect  ratio  10. 
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Grid 

ko  =  1,  M  = 

13 

k0  =  10,  M  = 

39 

k0  =  25,  M  = 

73 

||U''  -U2'1  H  co 

rate 

||uh  -  U2,*||co 

rate 

IK  -  U^Hoo 

rate 

64  x  64 

2.448706e+2 

— 

2.323499  e+4 

— 

2.151660e+7 

— 

128  x  128 

4.433401e+2 

-0.8564 

1.156875e+4 

1.0061 

3.572665e+6 

2.5904 

256  x  256 

6.193218e+2 

-0.4823 

2.889933e+2 

5.3231 

1.535344e+4 

7.8623 

512  x  512 

1.050433e-2 

15.8474 

1.251498e-l 

11.1732 

4.447023e-l 

15.0754 

1024  x  1024 

5.264725e-3 

0.9966 

1.706648e-2 

2.8744 

3.105638e-2 

3.8399 

2048  x  2048 

6.643254e-4 

2.9864 

1.993312e-3 

3.0979 

1.968995e-3 

3.9794 

4096  x  4096 

2.408868e-5 

4.7855 

1.040010e-4 

4.2605 

2.088858e-4 

3.2367 

Table  3.22:  Sound-soft  scattering  of  a  plane  wave  with  incidence  angle  50°  about  an  ellipse 
with  aspect  ratio  10. 


Grid 

II 

< 

II 

O 

13 

k0  =  10,  M  = 

35 

k0  =  25,  M  = 

61 

||uh-U2h||co 

rate 

Hl^-U^Hoo 

rate 

K  —  U2h\\oo 

rate 

64  x  64 

1.181451e-2 

— 

1.316527  e+1 

— 

1.879091e+2 

— 

128  x  128 

3.146854e-4 

5.2305 

5.405421e-l 

4.6062 

7.405706e+l 

1.3433 

256  x  256 

1.566558e-5 

4.3282 

4.087686e-3 

7.0470 

3.811382e-l 

7.6022 

512  x  512 

1.011771e-6 

3.9526 

2.520713e-4 

4.0194 

2.007000e-2 

4.2472 

1024  x  1024 

6.381159e-8 

3.9869 

1.567661e-5 

4.0071 

1.242483e-3 

4.0137 

2048  x  2048 

4.023256e-9 

3.9874 

9.784086e-7 

4.0020 

7.724207e-5 

4.0077 

4096  x  4096 

2.545110e-10 

3.9826 

6.104984e-8 

4.0024 

4.834473e-6 

3.9980 

Table  3.23:  Sound-hard  scattering  of  a  plane  wave  with  incidence  angle  50°  about  an 
ellipse  with  aspect  ratio  2. 


Grid 

II 

II 

O 

14 

k0  =  10,  M  = 

43 

k0  =  25,  M  = 

79 

K  —  U2'*  co 

rate 

||uh-U2h||oo 

rate 

||uh-U2h||co 

rate 

64  x  64 

5.859690e-2 

— 

4.558268 

— 

1.053224e+2 

— 

128  x  128 

9.846534e-3 

2.5731 

1.448027e-l 

4.9763 

4.461393e+l 

1.2392 

256  x  256 

1.884702e-4 

5.7072 

8.659281e-3 

4.0637 

5.844125e-l 

6.2544 

512  x  512 

9.615561e-6 

4.2928 

2.363689e-4 

5.1951 

2.368270e-2 

4.6251 

1024  x  1024 

4.412894e-7 

4.4456 

1.470104e-5 

4.0071 

1.454638e-3 

4.0251 

2048  x  2048 

2.845780e-8 

3.9548 

9.188393e-7 

4.0000 

9.027833e-5 

4.0101 

4096  x  4096 

1.589844e-9 

4.1619 

5.934903e-8 

3.9525 

5.631228e-6 

4.0029 

Table  3.24:  Sound-hard  scattering  of  a  plane  wave  with  incidence  angle  0°  about  an  ellipse 
with  aspect  ratio  3. 
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Grid 

ko  =  1,  M  = 

13 

k0  =  10,  M  = 

35 

k0  =  25,  M  = 

61 

_8 

IN 

1 

rate 

Uh  -  U2'1  oo 

rate 

|uh  -uMi|00 

rate 

64  x  64 

3.381436 

— 

5.287685 

— 

8.812860e+4 

— 

128  x  128 

1.076195 

1.6517 

1.602327 

1.7225 

6.882858e+l 

10.3224 

256  x  256 

1.593996e-l 

2.7552 

5.157602e-l 

1.6354 

2.620046 

4.7153 

512  x  512 

1.921666e-3 

6.3741 

9.974005e-3 

5.6924 

3.826752e-2 

6.0973 

1024  x  1024 

3.456720e-5 

5.7968 

2.426475e-4 

5.3612 

1.628882e-3 

4.5542 

2048  x  2048 

3.522082e-6 

3.2949 

1.769029e-5 

3.7778 

1.062220e-4 

3.9387 

4096  x  4096 

1.822888e-7 

4.2721 

9.543673e-7 

4.2123 

6.264534e-6 

4.0837 

Table  3.25:  Sound-hard  scattering  of  a  plane  wave  with  incidence  angle  50°  about  an 
ellipse  with  aspect  ratio  5. 


We  also  emphasize  that  the  scheme  converges  with  the  same  design  rate  for  all 
angles  of  incidence,  all  wavenumbers,  and  all  aspect  ratios.  The  actual  values  of  the 
error  may,  of  course,  depend  on  the  specific  parameters  involved.  For  example,  from 
Tables  3.16  through  3.25  one  can  see  that  as  the  wavenumber  ko  increases  while  all 
other  parameters  remain  the  same  (the  aspect  ratio,  the  grid,  etc.),  the  error  also 
increases  (maximum  norm  evaluated  across  the  domain).  On  the  other  hand,  the 
angle  of  incidence  does  not  affect  the  convergence  rate  and  does  not  noticeably  affect 
the  actual  error  either.  In  Figures  3.7,  3.8  and  3.9  we  show  the  dependence  of  the 
error  on  the  angle  of  incidence  for  both  sound-soft  (Dirichlet  boundary  condition) 
and  sound-hard  (Neumann  boundary  condition)  scattering  about  an  ellipse  of  aspect 
ratio  3.  We  see  that  for  both  ko  =  3,/co  =  15  and  ko  =  30  the  error  changes  by  less 
than  a  factor  of  2  over  the  entire  90°  range. 

For  the  same  setting  as  before,  we  also  conducted  a  series  of  computations  that 
corroborate  the  pollution  effect.  As  mentioned  in  Section  1.2,  to  maintain  the  same 
level  of  error  for  different  values  of  the  wavenumber  k ,  the  quantity  hnkn+1  must 
remain  constant,  where  h  is  the  grid  size  and  n  is  the  order  of  accuracy.  This 
means,  in  particular,  that  if  the  grid  size  h  is  halved,  then  the  same  level  of  error 
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X  1<f 


Figure  3.7:  Error  vs.  the  angle  of  incidence  for  sound-soft  (Dirichlet  BC)  and  sound-hard  (Neu¬ 
mann  BC)  scattering  about  an  ellipse  with  aspect  ratio  3  for  the  wavenumber  fco  =  3  and  the 
dimension  of  the  basis  (2.38)  M  =  21  computed  on  the  polar  grid  of  dimension  1024  x  1024. 

n 

shall  be  expected  if  the  wavenumber  k  increases  by  a  factor  of  2n+1 .  Specifically,  the 

2 

wavenumber  k  should  increase  by  a  factor  of  23  ps  1.5874  for  the  central  difference 

4 

second  order  scheme,  by  a  factor  of  2  s  ss  1.7411  for  the  fourth  order  compact 
scheme  (1.8),  and  by  a  factor  of  2e  «  1.7818  for  the  sixth  order  compact  scheme  of 
].  We  see  that  this  factor  is  always  less  than  two,  which  means  that  the  “points- 
per-wavelength”  quantity  does  not  stay  constant  but  rather  increases  as  k  grows. 
However,  the  higher  the  accuracy  n  of  the  scheme  the  slower  the  error  increases.  In 
Table  3.26,  we  summarize  the  computational  data  for  all  three  schemes.  Specifically, 
we  fix  the  value  of  k  =  4  for  the  grid  of  dimension  256  x  256,  and  then  vary  h  and 

n 

change  k  according  to  the  power  law  2 n+ 1  that  corresponds  to  each  scheme. 

We  next  perform  complexity  profiling  of  the  algorithm.  We  stress  that  our 
methodology  is  naturally  very  effective  for  problems  with  several  impinging  waves 
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Figure  3.8:  Error  vs.  the  angle  of  incidence  for  sound-soft  (Dirichlet  BC)  and  sound-hard  (Neu¬ 
mann  BC)  scattering  about  an  ellipse  with  aspect  ratio  3  for  the  wavenumber  ko  =  15  and  the 
dimension  of  the  basis  (2.38)  M  =  52  computed  on  the  polar  grid  of  dimension  1024  x  1024. 


x  10~3 


(a)  k0  =  30,  M  =  86 

Figure  3.9:  Error  vs.  the  angle  of  incidence  for  sound-soft  (Dirichlet  BC)  and  sound-hard  (Neu¬ 
mann  BC)  scattering  about  an  ellipse  with  aspect  ratio  3  for  the  wavenumber  ko  =  30  and  the 
dimension  of  the  basis  (2.38)  M  =  86  computed  on  the  polar  grid  of  dimension  1024  x  1024. 
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Grid 

k 

M 

Problem 

Diriclet 

Neumann 

||uh  —  u2,l||oo 

||u'‘-u2'l||oc, 

64  x  64 

1.3195 

14 

5.088513e-3 

1.107942e-3 

128  x  128 

2.2974 

18 

4.707891e-4 

1.385318e-4 

256  x  256 

4.0000 

24 

5.628050e-5 

1.061957e-5 

512  x  512 

6.9644 

32 

4.670625e-5 

4.257857e-5 

1024  x  1024 

12.125 

46 

4.203750e-5 

3.942691e-5 

2048  x  2048 

21.112 

66 

3.922412e-5 

3.755659e-5 

4096  x  4096 

36.758 

100 

3.764344e-5 

3.643574e-5 

Table  3.26:  Behavior  of  the  schemes  for  various  k  —  manifestation  of  the  pollution  effect. 


(see  discussion  in  one  dimensional  example  in  Chapter  2  and  also  Algorithm  11). 
This  is  useful  in  many  practical  problems,  therefore  information  about  complexity 
in  solving  a  scattering  problem  with  multiple  incident  angles  is  presented. 

Our  implementation  was  written  in  MATLAB,  and  the  linear  systems  obtained 
from  our  scheme  are  solved  via  MATLAB’s  built-in  direct  sparse  solver.  The  com¬ 
putations  were  performed  on  a  2.93  GHz  Quad-Core  Intel  Xeon  with  32  Gb  of  RAM 
running  on  Mac  OS  X. 

In  Table  3.27  we  present  sound-soft  scattering  of  a  plane  wave  about  an  ellipse 
with  an  aspect  ratio  of  2  for  the  wavenumber  k$  =  20.  The  dimension  of  the  basis 
(2.38)  is  M  =  43.  We  first  solve  for  an  incidence  angle  of  0°  and  then  use  the  obtained 
information  to  solve  the  same  problem  with  an  impinging  wave  at  a  different  incident 
angle  of  5°.  Column  4  shows  that  the  computational  complexity  of  the  method  is 
somewhat  better  than  linear  as  the  grid  is  refined.  This  behavior  is  caused  by  the 
direct  solver  used  which  was  developed  in  [  ].  Column  7  shows  the  linear  behavior 
of  the  computational  complexity  of  the  calculation  of  the  second  impinging  wave. 
The  matrices  Q  (2.64)  computed  for  the  first  impinging  wave  don’t  need  to  be 
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recalculated,  and  therefore  the  reduced  time  at  this  stage  is  due  to  the  calculation  of 
the  coefficients  and  obtaining  the  final  solution  by  solving  another  auxiliary  problem 
one  more  time  using  the  solver  from  [  >],  see  Section  2.3.6.  Column  8  shows  the  time 
saved  when  a  second  problem  is  solved  relative  to  the  first  one.  Similar  savings 
apply  for  any  number  of  different  incident  waves.  A  sound-hard  scattering  with  the 
same  setting  is  presented  in  Table  3.28.  We  find  a  similar  linear  behavior  and  a 
similar  saving  of  time  for  the  second  problem. 


Grid 

1  wave  6»(inc)  = 

0° 

2  wave  0(inc)  = 

5° 

Saving 

I|uh-u2h[|0O 

Time(s) 

Scaling 

||tU-U2'Too 

Time(s) 

Scaling 

64  x  64 

1.354291e+l 

6.2004e-l 

— 

1.110603e+l 

8.7971e-2 

— 

7.05 

128  x  128 

1.234442e-l 

1.4990 

2.4176 

1.317353e-l 

1.9846e-l 

2.2560 

7.55 

256  x  256 

6.971178e-3 

4.7169 

3.1466 

7.43461  le-3 

4.0498e-l 

2.0405 

11.65 

512  x  512 

4.245374e-4 

1.7626e+l 

3.7368 

4.521149e-4 

1.1563 

2.8553 

15.24 

1024  x  1024 

2.635393e-5 

6.6883e+l 

3.7944 

2.805860e-5 

3.9594 

3.4240 

16.89 

2048  x  2048 

1.644256e-6 

2.6894e+2 

4.0211 

1.750468e-6 

1.6856e+l 

4.2571 

15.96 

4096  x  4096 

1.027215e-7 

1.4513e+3 

5.3963 

1.093551e-7 

8.7533e+l 

5.1929 

16.58 

Table  3.27: 

CPU  times  for  sound-soft  scattering  of  a  plane  wave  with  incidence  angles  0° 

and  5°  about 
basis  M  =  43. 

an  ellipse  with  aspect  ratio  2  for  wavenumber  kt 

3  =  20.  The  dimension  of  the 

Grid 

1  wave  0(inc3  = 

0° 

2  wave  6>(inc)  = 

5° 

Saving 

||u'‘-u2h[|00 

Time(s) 

Scaling 

||tU-U2'd|oo 

Time(s) 

Scaling 

64  x  64 

6.366639 

6.1528e-l 

— 

5.899738 

1.3213e-l 

— 

4.66 

128  x  128 

5.268134  e-1 

1.4878 

2.4180 

1.096330 

2.0253e-l 

1.5328 

7.35 

256  x  256 

7.056579e-3 

4.7161 

3.1700 

7.205329e-3 

4.0841e-l 

2.0166 

11.55 

512  x  512 

4.305407e-4 

1.7395e+l 

3.6884 

4.414049e-4 

1.1562 

2.8309 

15.05 

1024  x  1024 

2.680923e-5 

6.6856e+l 

3.8434 

2.748845e-5 

3.9751 

3.4381 

16.82 

2048  x  2048 

1.673873e-6 

2.6953e+2 

4.0315 

1.715142e-6 

1.6876e+l 

4.2455 

15.97 

4096  x  4096 

1.046800e-7 

1.4521e+3 

5.3875 

1.071411e-7 

8.7569e+l 

5.1888 

16.58 

Table  3.28:  CPU  times  for  hard-soft  scattering  of  a  plane  wave  with  incidence  angles  0° 
and  5°  about  an  ellipse  with  aspect  ratio  2  for  wavenumber  fco  =  20.  The  dimension  of  the 
basis  M  =  43. 
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3.3  Transmission— Reflection  problems 

3.3.1  Piecewise  constant  coefficients 

The  numerical  simulation  of  the  simultaneous  transmission  and  scattering  of  waves 
off  a  given  shape  (an  ellipse)  is  done  using  a  computational  framework  similar  to  that 
of  Section  3.2,  except  that  instead  of  setting  a  boundary  condition  on  T  we  assume 
that  the  medium  inside  the  ellipse  is  characterized  by  a  constant  wavenumber  k\ 
(typically,  k\  >  ko).  At  the  interface  T  the  overall  solution  and  its  first  normal 
derivative  are  continuous. 

The  exterior  AP  and  its  discretization  remain  the  same  as  in  Section  3.2,  while 
the  interior  AP  is  formulated  on  the  rectangle  [—a  —  0.2,  a  +  0.2]  x  [— b  —  0.2,  b  +  0.2], 
where  a  and  b  are  the  major  and  minor  semi-axes  of  the  ellipse,  respectively.  We 
keep  a  =  1.8  and  vary  b  between  0.9  and  0.15,  which  yields  aspect  ratios  between  2 
and  12.  The  boundary  conditions  for  the  interior  AP  are  zero  Dirichlet  at  the  two 
horizontal  sides  of  the  rectangle,  and  local  Sonmierfeld-type  conditions  (complex)  at 
its  two  vertical  sides  (see  (2.34)).  The  latter  guarantees  the  unique  solvability  of  the 
interior  AP  (no  resonances),  see  [9,  Section  4.2]  or  [  ,  Section  5.2],  The  interior  AP 

is  discretized  by  the  compact  scheme  [  ]  with  fourth  order  accuracy  on  a  uniform, 
in  each  coordinate  direction,  Cartesian  grid.  It  is  then  solved  by  a  sparse  direct 
linear  solver  built  into  MATLAB.  To  simplify  the  monitoring  and  analysis  of  the 
grid  convergence,  the  grid  dimensions  for  the  interior  and  exterior  AP  are  always 
kept  the  same,  i.e. ,  those  two  grids  are  refined  synchronously.  As  in  Section  3.2,  the 
convergence  is  assessed  by  evaluating  the  maximum  volume  norm  of  the  difference 
between  the  numerical  solutions  obtained  on  two  consecutive  grids.  In  this  section 
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though,  it  is  done  independently  for  the  exterior  and  interior  parts  of  the  overall 
solution. 


Grid  Exterior 

k0  =  1,  fei  =  3,M  =  18  k0  =  5,  fci  =  15,  M  =  43  k0  =  10,  k0  =  30,  M  =  70 


u'‘u2h||oo  rate  ||uh  —  u2,l||oo  rate  ||uh  —  u2,l||oo  rate 


64  x  64 

4.595301e-3 

— 

3.193901 

— 

6.548224e+2 

— 

128  x  128 

1.680380e-4 

4.7733 

3.294341e-l 

3.2773 

1.939342 

8.3994 

256  x  256 

9.570738e-6 

4.1340 

1.892338e-3 

7.4437 

1.929997e-l 

3.3289 

512  x  512 

4.988526e-7 

4.2619 

1.133776e-4 

4.0610 

6.944723e-3 

4.7965 

1024  x  1024 

2.966397e-8 

4.0718 

6.982484e-6 

4.0213 

4.325129e-4 

4.0051 

2048  x  2048 

1.713626e-9 

4.1136 

4.349789e-7 

4.0047 

2.695934e-5 

4.0039 

Interior 

Ht^-U^Hoo 

rate 

||ufc  —  U2h||oo 

rate 

||uh-uah||oo 

rate 

64  x  64 

1.580527e-2 

— 

8.622188 

— 

7.361263e+2 

— 

128  x  128 

1.371018e-4 

6.8490 

4.000175e-l 

4.4299 

4.175908 

7.4617 

256  x  256 

8.490845e-6 

4.0132 

2.746052e-3 

7.1866 

2.013473e-l 

4.3743 

512  x  512 

4.202829e-7 

4.3365 

1.663763e-4 

4.0448 

8.816271e-3 

4.5134 

1024  x  1024 

2.442490e-8 

4.1049 

1.028977e-5 

4.0152 

5.509656e-4 

4.0001 

2048  x  2048 

1.366547e-9 

4.1597 

6.430838e-7 

4.0001 

3.435270e-5 

4.0035 

Table  3.29:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  aspect  ratio  2. 


The  data  Tables  3.29,  3.30  and  3.31  demonstrate  the  grid  convergence  for  two 
particular  sets  of  parameters  involved.  The  convergence  for  other  cases  that  we 
have  tried  with  piecewise  constant  k  looks  similar.  In  addition  to  showing  the 
convergence  data  in  Tables  3.29,  3.30  and  3.31,  we  also  plot  an  absolute  value  and 
real  and  imaginary  parts  of  the  solutions  that  we  have  computed,  see  Figures  3.10 
through  3.15. 
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Grid  Exterior 

fc0  =  l,fci  =  3,M  =  18  fc0  =  5,  fci  =  15,  M  =  43  fc0  =  10,  k0  =  30,  M  =  70 
||uh  —  u2h||oo  rate  |u,!'  —  u2h  |  |oo  rate  \uh  —  u2h||oo  rate 

64  x  64  3.397953e-2  —  1.928982  —  1.564833e+4 

128  x  128  5.257675e-4  6.0141  1.222531  0.6580  2.540071  12.5889 

256  x  256  2.968034e-5  4.1468  1.256795e-2  6.6040  2.505713  0.0196 

512  x  512  1.621693e-6  4.1939  7.215245e-4  4.1226  5.294069e-2  5.5647 

1024  x  1024  8.319173e-8  4.2849  4.686281e-5  3.9445  3.040675e-3  4.1219 

2048  x  2048  5.185747e-9  4.0038  2.980399e-6  3.9749  1.824703e-4  4.0587 


Interior 

||uh  _  u2h\\oo  rate  |u,!'  —  u2h | |oo  rate  \uh  —  u2h||oo  rate 

64  x  64  5.915841e-2  —  4.502945  —  2.997456e+4 

128  x  128  4.727249e-4  6.9674  1.476710  1.6085  1.454990e+l  11.0085 

256  x  256  9.910096e-6  5.5760  7.240317e-2  4.3502  3.219025  2.1763 

512  x  512  8.298524e-7  3.5780  7.604703e-4  6.5730  5.533717e-2  5.8622 

1024  x  1024  3.405935e-8  4.6067  5.105078e-5  3.8969  3.352915e-3  4.0448 


2048  x  2048  2.025545e-9  4.0717  3.330358e-6  3.9382  2.093157e-4  4.0017 


Table  3.30:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  aspect  ratio  3. 


Grid  Exterior 

k0  =  1,  fci  =  3,  M  =  17  fc0  =  5,  fci  =  15,  M  =  42  k0  =  10,  k0  =  30,  M  =  68 
||uh  —  u2h||oo  rate  ||uh  —  u2,l||oo  rate  |uh  —  u2h\\oo  rate 

128  x  128  3.339419  —  1.011317e+2  —  8.966609e+2 

256  x  256  2.525248e-3  10.3690  4.429088  4.5131  1.265494e+l  6.1468 

512  x  512  3.016655e-4  3.0654  3.727909e-2  6.8925  5.603543e-l  4.4972 

1024  x  1024  7.296313e-5  2.0477  4.070004e-3  3.1953  6.315946e-3  6.4712 

2048  x  2048  9.046537e-6  3.0117  3.237602e-4  3.6520  3.668013e-4  4.1059 


Interior 

||uh  —  u2h||oo  rate  |uh  —  u2,l||oo  rate  ||uh  —  u2,I||00  rate 

128  x  128  2.895946  —  7.446942e+2  —  4.085428e+4 

256  x  256  2.649119e-3  10.0943  4.907694  7.2455  1.662587e+2  7.9409 

512  x  512  2.590272e-4  3.3543  4.386424e-2  6.8059  1.008079  7.3657 

1024  x  1024  5.793892e-5  2.1605  3.933254e-3  3.4792  8.501939e-3  6.8896 


2048  x  2048  6.983459e-6  3.0525  3.009321e-4  3.7082  3.740583e-4  4.5065 


Table  3.31:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  aspect  ratio  12. 
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total  field  abs 


(a)  Aspect  ratio  2 

total  field  abs 


(b)  Aspect  ratio  10 


Figure  3.10:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  k\  =  20  (inside)  and  k0  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  1024  x  1024  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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(a)  Aspect  ratio  2 


total  lietd.  teal 


(b)  Aspect  ratio  10 


Figure  3.11:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  k\  =  20  (inside)  and  fc0  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  1024  x  1024  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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(a)  Aspect  ratio  2 

total  field  imag 
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(b)  Aspect  ratio  10 


Figure  3.12:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about 
an  ellipse  with  k\  =  20  (inside)  and  fc0  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  1024  x  1024  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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total  Held,  abs 


(a)  Aspect  ratio  3 

total  field,  abs 


(b)  Aspect  ratio  5 


Figure  3.13:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  180°  about 
an  ellipse  with  k\  =  20  (inside)  and  fc0  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  513  x  513  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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total  field,  real 


(a)  Aspect  ratio  3 

total  field,  real 


(b)  Aspect  ratio  5 


Figure  3.14:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  180°  about 
an  ellipse  with  k\  =  20  (inside)  and  ko  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  513  x  513  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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total  field,  imag 


(b)  Aspect  ratio  5 


Figure  3.15:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  180°  about 
an  ellipse  with  k\  =  20  (inside)  and  fc0  =  10  (outside).  Absolute  value  of  the  total  field  is 
shown.  The  grid  dimension  is  513  x  513  for  both  the  interior  AP  (Cartesian)  and  exterior 
AP  (polar). 
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3.3.2  Piecewise  smooth  coefficients 


In  this  section,  we  keep  the  computational  setting  the  same  as  in  Section  3.3.1, 
except  that  we  allow  for  a  smooth  variation  of  the  wavenumber  inside  the  ellipse: 


k 


kie-W(r-r0)N^  if  r  ^  r0, 
k\,  if  r  >  r0, 


(3.8) 


where  r  =  ^/x2  +  y2  and  ro  =  1.6.  The  variable  coefficient  Helmholtz  equation 
is  approximated  with  fourth  order  accuracy  by  the  compact  scheme  of  [  ].  In  Ta¬ 
ble  3.32,  we  show  the  results  for  the  ellipse  with  aspect  ratio  2:  a  =  1.8  and  b  =  0.9. 
In  Table  3.33,  we  show  the  results  for  the  ellipse  with  aspect  ratio  3:  a  =  1.8  and 
b  =  0.6.  In  Table  3.34,  we  show  the  results  for  the  ellipse  with  aspect  ratio  5:  a  =  1.8 
and  b  =  0.36.  These  results  corroborate  the  design  fourth  order  convergence  rate  of 
the  algorithm  in  the  case  of  variable  coefficients. 


Grid  Exterior 

k0  =  1,  fei  =  3,M  =  50  k0  =  5,fci  =  15,  M  =  50  k0  =  10,  k0  =  30,  M  =  69 


u'l-u2h||00  rate  ||uh  —  u2,l||oo  rate  ||uh  —  u2,l||oo  rate 


128  x  128 

1.196903e-l 

— 

2.740357 

— 

9.679627 

— 

256  x  256 

8.888395e-4 

7.0732 

1.071586e-l 

4.6765 

8.093132e-l 

3.5802 

512  x  512 

5.054118e-5 

4.1364 

1.339804e-3 

6.3216 

8.933817e-3 

6.5013 

1024  x  1024 

9.918967e-7 

5.6711 

4.031620e-5 

5.0545 

5.392737e-4 

4.0502 

2048  x  2048 

7.239554e-8 

3.7762 

2.729896e-6 

3.8844 

4.167435e-5 

3.6938 

Interior 

||ufc  —  Uah||oo 

rate 

||uh-U2,l||oo 

rate 

l|uh-u2'l||00 

rate 

128  x  128 

2.169952e-l 

— 

5.413927 

— 

2.473694e+l 

— 

256  x  256 

1.267208e-3 

7.4199 

1.422706e-l 

5.2500 

1.466237 

4.0765 

512  x  512 

4.573756e-5 

4.7921 

1.506757e-3 

6.5610 

5.642709e-2 

4.6996 

1024  x  1024 

6.765407e-7 

6.0791 

2.968551e-5 

5.6655 

4.756129e-4 

6.8905 

2048  x  2048 

5.114153e-8 

3.7256 

2.285437e-6 

3.6992 

2.445819e-5 

4.2814 

Table  3.32:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about  an 
inhomogeneous  ellipse  with  aspect  ratio  2  and  interior  wavenumber  given  by  formula  (  3.8). 
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Grid  Exterior 

k0  =  l,fcr  as  3,M  =  44  k0  =  5,  fci  =  15,  M  =  4.3  k0  =  10,  k0  =  30,  M  =  70 


—  u2^1  1 1  oo  rate  ||uh  —  u2h||oo  rate  ||uh  —  u2h||oo  rate 


128  x  128 

1.342633e-3 

— 

3.189426e-l 

— 

2.273503 

— 

256  x  256 

4.704559e-5 

4.8349 

2.112324e-3 

7.2383 

1.428778 

0.6701 

512  x  512 

1.972406e-6 

4.5760 

1.037815e-4 

4.3472 

8.515279e-3 

7.3905 

1024  x  1024 

1.128138e-7 

4.1279 

6.701632e-6 

3.9529 

4.664540e-4 

4.1902 

2048  x  2048 

6.620319e-9 

4.0909 

5.703865e-7 

3.5545 

2.867676e-5 

4.0238 

Interior 

||uh-uah||oo 

rate 

||ufc  —  Uafc||oo 

rate 

Hu^-U^lloc 

rate 

128  x  128 

3.067343e-3 

— 

1.043931 

— 

9.253122 

— 

256  x  256 

6.693782e-5 

5.5180 

3.225823e-2 

5.0162 

2.960923 

1.6439 

512  x  512 

1.396376e-6 

5.5831 

1.002967e-4 

8.3292 

4.859728e-2 

5.9290 

1024  x  1024 

8.178574e-8 

4.0937 

6.580268e-6 

3.9300 

4.749737e-4 

6.6769 

2048  x  2048 

5.333208e-9 

3.9388 

4.146442e-7 

3.9882 

2.935093e-5 

4.0164 

Table  3.33:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about  an 
inhomogeneous  ellipse  with  aspect  ratio  3  and  interior  wavenumber  given  by  formula  (  3.8). 


Grid  Exterior 

fco  =  l,fei  =  3,M  =  50  fc0  =  5,fei  =  15,M  =  50  k0  =  10,  k0  =  30,  M  =  69 


U^-u^Hoo  rate  ||uh  —  u2,l||oo  rate  | \uh  —  u2,l||oo  rate 


128  x  128 

1.196903e-l 

— 

2.740357 

— 

9.679627 

— 

256  x  256 

8.888395e-4 

7.0732 

1.071586e-l 

4.6765 

8.093132e-l 

3.5802 

512  x  512 

5.054118e-5 

4.1364 

1.339804e-3 

6.3216 

8.933817e-3 

6.5013 

1024  x  1024 

9.918967e-7 

5.6711 

4.031620e-5 

5.0545 

5.392737e-4 

4.0502 

2048  x  2048 

7.239554e-8 

3.7762 

2.729896e-6 

3.8844 

4.167435e-5 

3.6938 

Interior 

||uh-u2'i||00 

rate 

Uh  -  U2'*  oo 

rate 

||uh  —  U2  ^  1 1  oo 

rate 

128  x  128 

2.169952e-l 

— 

5.413927 

— 

2.473694e+l 

— 

256  x  256 

1.267208e-3 

7.4199 

1.422706e-l 

5.2500 

1.466237 

4.0765 

512  x  512 

4.573756e-5 

4.7921 

1.506757e-3 

6.5610 

5.642709e-2 

4.6996 

1024  x  1024 

6.765407e-7 

6.0791 

2.968551e-5 

5.6655 

4.756129e-4 

6.8905 

2048  x  2048 

5.114153e-8 

3.7256 

2.285437e-6 

3.6992 

2.445819e-5 

4.2814 

Table  3.34:  Transmission  and  scattering  of  a  plane  wave  with  incidence  angle  40°  about  an 
inhomogeneous  ellipse  with  aspect  ratio  5  and  interior  wavenumber  given  by  formula  (  3.8). 
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Chapter  4 

Discussion  and  Conclusion 


4.1  Discussion 

We  have  described  a  combined  implementation  of  the  method  of  difference  poten¬ 
tials  along  with  the  compact  high  order  accurate  finite  difference  schemes  for  the 
numerical  solution  of  wave  propagation  problems  in  the  frequency  domain.  The 
governing  Helmholtz  equation  is  approximated  on  a  regular  structured  grid,  which 
is  efficient  and  entails  a  low  computational  complexity.  At  the  same  time,  the 
method  guarantees  no  loss  of  accuracy  for  curvilinear  non-conforming  boundaries, 
and  can  also  handle  variable  coefficients  that  describe  a  non-homogeneous  medium. 
As  such,  the  resulting  method  provides  a  viable  alternative  to  both  BEM  and  high 
order  FEM. 

The  performance  of  the  method  and,  in  particular,  its  design  high  order  accu¬ 
racy,  has  been  corroborated  numerically  by  solving  a  variety  of  2D  interior,  exterior, 
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and  reflection/transmission  Helmholtz  problems,  including  those  with  variable  co¬ 
efficients,  on  Cartesian  and  Polar  grids  for  non-conforming  boundaries/interfaces 
shaped  as  circles  and  ellipses. 

Among  other  advantages  of  the  proposed  methodology  is  its  capability  to  accu¬ 
rately  reconstruct  the  solution  and/or  its  normal  derivative  directly  at  the  boundary. 
This  is  done  without  interpolating  and/or  using  one-sided  differences,  such  as  in  con¬ 
ventional  FD,  and  with  no  additional  developments  needed  as  in  FEM,  see,  e.g.,  [  ]. 

Additional  advantages  of  the  method  are  the  absence  of  any  singular  integrals  or 
similar  constructs,  the  minimum  number  of  unknowns  that  characterize  the  discrete 
solution  —  just  one  per  grid  node,  and  the  same  number  of  boundary  conditions 
needed  for  the  scheme  as  that  needed  for  the  underlying  differential  equation. 

For  exterior  problems  we  constructed  auxiliary  problems  (see  Section  2. 3. 2. 2) 
with  the  appropriate  artificial  boundary  conditions  (ABCs),  see  [  ].  For  the  con¬ 

stant  coefficient  2D  Helmholtz  equation  (typical  for  the  far  field),  the  APs  were 
formulated  using  polar  coordinates,  which  enables  a  natural  and  efficient  implemen¬ 
tation  of  the  exact  non-local  ABCs  in  Fourier  space.  The  proposed  methodology  is 
extended  for  the  combined  reflection/transmission  problems.  The  latter  formulation 
involves  a  joint  solution  of  the  interior  and  exterior  Calderon’s  BEPs  constructed  at 
the  interface  between  the  interior  and  exterior  sub- regions. 

Thus  far,  we  have  computed  solutions  only  for  circular  and  elliptical  boundaries. 
The  case  of  general  smooth  boundaries  was  analyzed  theoretically  in  [  12,  Appendix 
A],  and  the  corresponding  Taylor-based  extension  operators  have  been  developed. 
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4.2  Conclusions  and  future  research 

We  have  simulated  a  broad  range  of  constant  and  variable  coefficient  2D  test  cases 
for  non-conforming  boundaries/interfaces  on  regular  structured  grids.  Our  com¬ 
putations  convincingly  corroborate  the  design  high  order  accuracy  of  the  proposed 
method. 

A  study  is  needed  to  explore  alternative  strategies  for  choosing  M  -  the  dimen¬ 
sion  of  the  basis  used  for  representing  the  solution  at  the  boundary  T  as  in  (2.38)  or 
(2.59).  In  addition  to  that,  the  possibility  of  using  other  bases  has  to  be  examined. 
The  ideas  of  reduced  order  modeling  [  ]  (or,  similarly,  principal  component  analysis 

or  proper  orthogonal  decomposition)  has  to  be  investigated  to  further  reduce  the 
dimension  of  the  basis  on  T.  Other  bases  may  be  more  convenient  to  use  like  in  the 
case  of  piece- wise  parametrization  [  ]. 

Actual  computations  of  the  scattering  and  transmission/scattering  solutions  for 
arbitrarily  shaped  domains  (other  than  circles  or  ellipses)  will  need  to  be  performed. 
The  corresponding  general  construct  of  the  equation  based  extension  has  been  de¬ 
veloped  in  [42,  Appendix  A], 

The  domains  of  a  more  general  shape  (beyond  circles  and  ellipses)  have  to  be 
developed.  Arbitrary  smooth  boundaries  will  require  a  more  general  construction 
of  the  extension  operators  (see  Appendix  in  [  ]). 

Yet  another  direction  for  future  work  will  be  to  allow  for  multiple  sub-regions, 
for  example,  multiple  scatterers  immersed  into  the  same  background  medium.  The 
simplest  case  will  amount  to  considering  a  piece-wise  constant  function  k2(x,y)  in 
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the  Helmholtz  equation  (1.1),  while  more  elaborate  settings  may  also  include  the 
variation  of  material  characteristics. 

The  proposed  methodology  needs  to  be  extended  to  3D  wave  propagation  which 
requires  that  the  coordinates  associated  with  a  curve  be  replaced  with  surface- 
oriented  coordinates,  see  [  ]. 

Finally,  the  extension  of  the  proposed  methodology  to  time-dependent  prob¬ 
lems  (e.g.,  the  wave,  i.e.,  d’Alembert,  equation  instead  of  the  Helmholtz  equation) 
requires  additional  theoretical  developments. 

Note,  that  even  though  the  current  implementation  and  discussion  focuses  on 
wave  propagation  problems,  the  method  of  difference  potentials  is  capable  of  ad¬ 
dressing  a  considerably  broader  range  of  formulations,  including  problems  in  heat 
transfer,  elasticity,  fluid  dynamics,  and  other  areas,  see  for  example  [  ]. 
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